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ABSTRACT 

The observed high covering fractions of neutral hydrogen (Hi) with column densi- 
ties above ~ 10^^ cm“^ around Lyman-Break Galaxies (LBGs) and bright quasars 
at redshifts z ^ 2 — 3 has been identified as a challenge for simulations of galaxy 
formation. We use the EAGLE cosmological, hydrodynamical simulation, which has 
been shown to reproduce a wide range of galaxy properties and for which the subgrid 
feedback was calibrated without considering gas properties, to study the distribution 
of Hi around high-redshift galaxies. We predict the covering fractions of strong Hi 
absorbers (A^hi ^ lO^^cm”^) inside haloes to increase rapidly with redshift but to 
depend only weakly on halo mass. For massive (M 200 ^ 10^^ Mq) halos the cover¬ 
ing fraction profiles are nearly scale-invariant and we provide fitting functions that 
reproduce the simulation results. While efficient feedback is required to increase the 
Hi covering fractions to the high observed values, the distribution of strong absorbers 
in and around halos of a fixed mass is insensitive to factor of two variations in the 
strength of the stellar feedback. In contrast, at fixed stellar mass the predicted Hi dis¬ 
tribution is highly sensitive to the feedback efficiency. The fiducial EAGLE simulation 
reproduces both the observed global column density distribution function of Hi and 
the observed radial covering fraction profiles of strong Hi absorbers around LBGs and 
bright quasars. 

Key words: methods: numerical - galaxies: formation - galaxies: high-redshift - 
galaxies: absorption lines - quasars: absorption lines - intergalactic medium 


1 INTRODUCTION 

Galaxies need to acquire large quantities of fresh gas 
from the intergalactic medium (IGM) t o sustain their star- 
form ation activities through time (e.g., iBauermeister et al.l 
I2OIOI) . Simulations predict that a large fraction of the 
accreting material enters halos relatively cold and could 
therefore contain significant amounts of neutr al gas (e.g., 
iFumagalli et al.l I2OIII : Ivan de Voort et al.l I2OI2I ) . The pres¬ 
ence of shock-heated gas complicates the journey of the ac¬ 
creting gas onto galaxies. Moreover, energetic feedback from 
stars and active galactic nuclei (AGN), which regulate the 
consumption of the accreted gas and launch galactic out¬ 
flows, affect the dynamics and chemical composition of gas 
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around galaxies. As a result, the complex distribution of gas 
around galaxies contains the finger-prints of the aforemen¬ 
tioned processes and studying it, mostly by analysing the 
absorption signature of neutral hydrogen and metals in the 
spectra of bright background sources, is of great value for 
understanding galaxies and the physical processes that reg¬ 
ulate them. 


Observations and simulations show that both the num¬ 
ber of Hi absorbers and their t ypical column densities 
increase closer to galax i es (e.g.. 


_ _ — - __ lAdelberger et al.l 20031: 

Chen&^dulcha^^ 200 ^ iRajdc et al.l I2OI2I : iTurneretall 


2014 Rahmati fc Scha!^ 20141 ). This suggests that ab¬ 


sorbers with higher Hi column densities are better probes 
of the gas in the vicinity of galaxies. Cosmological simula¬ 
tions, however, suggest that most strong Hi absorbers, such 
as Lyman Limit Systems (LLSs; with Ahi > 10^^'^ cm“^) 
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and Damped Lyman-a systems (DLAs; with A^hi ^ 
10^°'®cm“^), are close to galaxies th at are too faint to 
be ea sily detectable in current surveys l|Rahmati fc Schavel 
l2014h . which is in agreement with the lack of detected 
counterparts close to most of strong Hi absorbers (e.g., 
iFumaealli et al.ll2015l l. Not knowing the properties of the 
host galaxy complicates the use of strong Hi absorbers to 
study the relation between galaxies and their environments. 
This problem can, however, be circumvented by studying the 
distribution of Hi absorbers around easily detectable bright 
galaxies. 

Several modern observational campaigns have adopted 
this galaxy-centred approach by using quasar absorp- 


tion lines to systematically 
of neutral hvdrogen around 

investigate the distribution 
massive galaxies at different 

epochs fe.g.. Adelberger et al.l 20031: Hennawi et al.l 20061: 

Chen & Mulchaevll2009|: iRak 

c et al.l 2 OI 2 I: iTumlinson et al.l 

2 OI 3 I: Prochaska et al.l l2013a 

: Turner et al.l 2014 Il For in- 


stance, 

of Hi around Lyman-Break galaxies at z ~ 2 and found that 
there is a ~ 30% chance for finding LLSs in the spectra of 
background quasars that have impact parameters less than 
~ 100 proper kpc (hereafter pkpc). Noting that this impact 
parameter is comparable to the virial radii of LBGs at z ~ 2, 
this result implies that LLS have a covering fraction of 30% 
within the virial radius of LBGs. IProchaska et al.l (l2013bl l 
showed that the abundance of Hi absorbers is significantly 
enhanced out to several virial radii from bright quasars at 
a ~ 2. They found that there is a more than 60% chance 
of finding a LLS within 150 pkpc from a bright quasar 
dProchaska et al]l2013al L which is comparable to the typical 
virial radius of the halos with M 200 ~ 10^^ ® M© that are 
expected to host the observed quasars at z ^ 2. 

Motivated by recent observational constraints, sev- 


of Hi around galaxies (e. 

g., IFaucher-Giguere & Keres 

2 OIII: IFaucher-Giguere et al.l 

2 OI 5 I: IFumagalli et al.l 

2011, 

2 OI 4 I: IShen et al.l 

2 OI 3 I: Rakic et al.l 20131: Erkall 

2014; 

Meiksin et al.l 2014 

). However, reproducing the relatively 


large observed covering fractions of strong Hi absorption 
arou nd massive galaxies tu r ned out to be a major challe nge 
fe.g.. [Rimagalli et al.|[20l3 : lFaucher-Giguere et al.ll2015h . 

Previous studies of high column density Hi around mas¬ 
sive galaxies were based on the a nalysis of simulations that 
zoom into only one gal axy (e.g., iFaucher-Giguere fc Ker^ 
I 2 OIII: IShen et all |2013|) . or a handful of galaxies span- 
ning a limited range in mass and redshift (e.g ., 
iFaucher-Giguere et al.l T 2 OI 5 I : iFumagalli et al.l 1201 ll . 1201411 . 
Given the diversity of observed galaxies, one may expect 
a large intrinsic variation in the spatial distribution of Hi 
from one galaxy to another. Consequently, a large sample 
of simulated galaxies is required to predict the average dis¬ 
tribution of Hi and to compare it robustly with observa¬ 
tions. Because observational constraints are at present lim¬ 
ited to galaxies residing in relatively massive halos (M 200 ^ 
10^^ M 0 ), which are rare, especially at high redshifts, sim¬ 
ulating a large number of them requires large cosmolog¬ 
ical volumes. Moreover, without cosmological simulations 
it is not straightforward to check whether the simulations 
satisfy other important constraints on the cosmic distri¬ 
bution of Hi, for instance the global Hi column density 
distribution function and/or the Hi cosmic density, which 


has been reprod uced successfully i n recent cosmo l ogical 
simulations (e.g., Mt^^^aiJ_ 20n|;_ Mc^uirin_et_^ 20n • 


RnLmati et al. I l2013al : ID ave et al. I I2013I : IVogelsbergCT et al.l 

2014h . The aforementioned issues may limit the power of 


studies that use small numbers of zoom simulations and indi¬ 
cate that cosmological simulations of representative volumes 
are needed to study the average distribution of Hi around 
large numbers of galaxies. 

The strength and implementation of feedback mecha¬ 
nisms are also crucial factors for si mulations of the distri¬ 
bution of gas around g alaxies (e.g., IFaucher-Giguere et al.l 
l2015l : [Suresh et ^l2015l L For instance, galactic winds driven 
by stellar feedback can change the distribution of Hi around 
galaxies by carrying the cold neutral gas farther away from 
galaxies and by providing resistance against the accretion 
of gas. While feedback implementations vary widely, simu¬ 
lations that use strong stellar feedback have been more suc¬ 
cessful in reproduc ing the LLS covering fractions observed 
aroun d LBGs ie.g.. lShen et al.ir2013l : IFaucher-Giguere et al.l 
l2015ll . which are significant ly underproduced in simulations 
with weak feedback ( e.g., Faucher-Giguere fc Kerej I 2 OIII : 
IFumagalli et al.ll201ll . 12014 1. Moreover, feedback from ac¬ 
tive galactic nuclei (AGN), which is required to form rea¬ 
sonable galaxies in very massive haloes and is missing from 
most previous studies, may affect the observed Hi distri¬ 
bution in and around halos expected to host quasars at 
z ^ 2. However, several simulations indicate that only the 
strongest absorbers, which on average reside closer to or even 
inside galaxies, are s ignificantly affected by feedback (e.g.. 


Theuris_et_aL 20021: Altav et ^ I 2 OI 3 I : iRahmati fc Schavel 
2 OI 4 I : Bird et al. 2014). Hence, the more important effect 


may be that feedback changes the relation between stellar 
mass and halo mass and t hus the predicted Hi distribution 
at fixed stellar mass fe.g.. lRakic et al]|2013ll . 

In this work, we study the Hi distribution around galax¬ 
ies using state-of-the-art cosmological hydrodynamical simu¬ 
lations. For this purpose, we use the Evolution and Assembly 
of Galaxies and th eir Environment (EAGLE) simulations 
dSchave et akllioisl . hereafter S15). The large cosmological 
volume of the main EAGLE run (100 cMpc) together with 
its relatively high resolution for a simulation of this type, al¬ 
low us to study large numbers of halos with masses similar 
to those targeted by recent observations, without compro¬ 
mising the resolution needed to simulate the distribution of 
relevant Hi systems (e.g., LLSs). Efficient stellar and AGN 
feedback enables the simulation to successfully reproduce 
a large number of basic observed characteri stics of galax¬ 
ies over wide mass and redshift ranges (S1 5; iFurlong et al.l 
l2014l : [Schaller et al.ll201 j : ICrain et al.ll2015l l and S15 already 
showed that EAGLE also reproduces the observed present- 
day column density distributions of CIV and OVI. These 
factors make EAGLE ideal for studying the gas distribution 
around galaxies. 

We combine the EAGLE simu lations with the accurat e 
photoionization corrections from IRahmati et al.l (l2013al L 
which are based on high-resolution radiative transfer calcu¬ 
lations. After showing the success of the simulation in repro¬ 
ducing the observed cosmic distribution of Hi, we look at the 
Hi distribution around galaxies from z = 4 to z = 1, bracket¬ 
ing the era during which the cosmic star formation density 
peaked. We focus on the strong Hi absorbers whose high 
covering fractions were found to be difficult to reproduce by 
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previous simulations (e.R., Faucher-Giguer e & Keres 20lJ; 


iFaucher-Giguere et alj [^ISl : iFumagalli et aLn2011 . 2014ll . 
We predict that strong Hi absorbers, such as LLSs and 
DLAs, have a mean covering fraction within the virial ra¬ 
dius that increases rapidly with redshift, but depends only 
weakly on the halo mass (or star formation rate) at fixed 
redshift, suggesting that the distribution of absorbing gas 
around galaxies has a similar shape for different halo masses. 
Indeed, we show that the covering fraction of LLSs, sub 
DLAs and DLAs around massive galaxies (M 200 ^ 10^^ M©) 
follows profiles with similar shapes but different scale lengths 
that are tied to the virial radius and redshift. 

We construct samples of simulated galaxies by match¬ 
ing the halo masses and redshif t s to th ose used in the obser¬ 
vationa l studies of iRudie et al.l (l2012l l and IProchaska et al.l 
ll2013btl for LBGs and quasars, respectively. Accounting for 
the uncertainties in the amplitude of the Ultra Violet Back¬ 
ground (UVB) photoionization rate, our predictions are in 
excellent agreement with the observed Hi distributions. This 
shows that cosmological hydrodynamical simulations that 
are successful in reproducing reasonable galaxy properties, 
are also capable of predicting gas distributions in agree¬ 
ment with current observations. We conclude that there is 
no obvious missing ingredient in our general understanding 
of galaxy formation and evolution required to explain the 
observed Hi distributions around LBGs and bright quasars 
at z 2. 


The structure of this paper is as follows. In ^we intro¬ 
duce our cosmological simulations and discuss the photoion¬ 
ization corrections required for obtaining the Hi column den¬ 
sities and calculating their distribution around galaxies. We 
present our predictions for the Hi covering fractions and 
their evolution in m We compare the predictions with re¬ 
cent observations in Sjl] and discuss the impact of feedback 
on our results in lj5] We conclude in ^ 


2 SIMULATION TECHNIQUES 

In this section we briefly describe the hydrodynamical sim¬ 
ulations that we use to predict the Hi distributions. We fur¬ 
ther explain our halo finding method (am and our pho¬ 
toionization correction required for the Hi column density 
calculations (am. 


2.1 Hydrodynamical simulations 

We use the reference simulation of the EAGLE project, de¬ 
scribed in S15, as our fiducial simulation. The cosmolog¬ 
ical simulation was performed using a significantly modi¬ 
fied and extended version of the smoothed particle hydrody¬ 
nami cs (SPH) code GADGET-3 (last described in ISpringell 
l2005h . In particular, we use “Anarchy” (Dalla Vecchia, in 
preparation; see also appendix A of S15) which is an up¬ 
dated hydrodynamics algorithm incorp orating the pre ssure 
entropy formulation of SPH derived by Ho^kina 12013) and 
the time-step limiter of iDurier fc Dalla Vecchia ll2012tl (see 
Appendix n where the impact of using “Anarchy” is dis¬ 
cussed). The subgrid physics u sed in the simulatio n is based 
on that of the OWLS project (ISchave et akllioioh with nu¬ 
merous important improvements. Stellar and AGN feedback 
are implemented using the stochastic, thermal prescription 
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of lPalla Vecchia fc Scha(^ ll2012l f. without turning off radia¬ 
tive cooling or the hydrodynamics. Galactic winds develop 
naturally, without predetermined mass loading factors or ve¬ 
locities. We use a metallicity-dependent subgrid model for 
star formation together w ith the pressure-depend e nt sta r 
formation prescription of ISchave fc Dalla Vecchial ll2008ll . 
The feedback from AGN is updated such that the subgrid 
model for accretio n of gas onto black holes a ccounts for an¬ 
gular momentum iRosas- Guevara et al.l[201^ L A metallicity 
and density dependent stellar feedback efficiency is adopted 
to account, respectively, for greater thermal losses when the 
metallicity increases and for residu al spurious resolution de¬ 
pendent numerical radia tive losses llPalla Vecchia fc Schava 
I 2 OI 2 I : [Crain et aI]|2015l L The implementation of metal en¬ 
richment is similar to that of th e OWLS project and is de¬ 
scribed in IWiersma et af] ll2009bll. We fol l ow th e abundances 
of eleven elements assuming a Chabrierl (l2003tl initial mass 
function. These abundances are used for calculating radia¬ 
tive cooling/heating rates, element-by-element and in the 
pre sence of the uniform cosm ic microwave background and 
the iHaardt fc Madaul (1200 il l UVB model (IWiersma et al.l 
l2009al L The simulation is calibrated based on the present 
day observed galaxy stellar mass function and galaxy sizes, 
which are reproduced with unp recedented accura cy for a hy¬ 
drodynamical simulation ('S15: lGrain et al1l2015l L The same 
simulation also shows very good agreement with other ob¬ 
served galaxy properties such as the observed galaxy specific 
star formation rates, passive fractions, Tully-Fisher relation 
and the distribution of metals in the intergalactic medium 


Ischaller et al.l 

2 OI 4 IL the evolution of the galaxy stellar 

mass function ( 

Furlong et al.l 20141') and the molecular hv- 

droeen content of ealaxies (Lag;os et al.||2015fl. 


The adopted cosmological parameters are based on 
the most recent Planck results: {flm = 0.307, fib = 
0.04825, flA = 0.693, as = 0.828 8, = 0.9611, h = 

0.6777} dPlanck GollaborationI [2^013h . Our reference simu¬ 
lation, Ref-Ll DON 1504 , has a periodic box of L = 100 co¬ 
moving Mpc (cMpc) and contains 1504® dark matter par¬ 
ticles with mass 9.7 x 10® M© and an equal number of 
baryonic particles with initial mass 1.81 x 10® M©. The 
Plummer equivalent gravitational softening length is set to 
ecom = 2.66 comoving kpc (ckpc) and is limited to a maxi¬ 
mum physical scale of Cprop = 0.7 pkpc. We use simulations 
with different feedback implementations to test the impact 
of feedback variations on our results in ^ Simulations with 
different box sizes and resolutions are used to study the im¬ 
pact of those factors on our results in Appendix m Table m 
summarizes the simulations we used in this work. 


2.2 Identifying galaxies in EAGLE 

We identify galaxies using the Friends-of-Friends (FoF) al¬ 
gorithm to select groups of dark matter particles that are 
near each other (i.e., FoF halos), choosing a linking length 
of 6 = 0.2. In other words, we assume that galaxies reside in 
dark matter halos. In the next step, we group gravitationally 
bound particles of unique structures (subhalo s) using SUB¬ 
FIND dSpringel et al.l[200ll : iDolag et af]|2009ll . We identify 
the centre of each halo/galaxy as the position of the particle 
with the minimum gravitational potential in that halo. Then 
we define the virial radius, r 2 oo, as the radius within which 
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Table 1. List of cosm ological simulations used in this work. The first four simulations use model ingredients identical to the EAGLE 
reference simulation of ISchave et al.l ll2nifji . while the higher-resolution Recal-L025N0752 has been re-calibrated to the observed present- 
day galaxy mass function. Model NoAGN does not include AGN feedback, WeakFB and StrongFB use half and twice as strong stellar 
feedback compared to the reference simulation, respectively (see lCrain et al]|2015h . Models NoFB and NoFB-Gadget, do not include 
any feedback from stars or AGN. While the former uses “Anarchy” (Dalla Vecchia in prep) for the hydrodynamical calculations, the 
latter uses the standard GADGET-3 implementation jSDringe]||2005l l. as does model Ref-Gadget. From left to right the columns show: 
simulation identifier; comoving box size; number of particles (there are equally many baryonic and dark matter particles); initial baryonic 
particle mass; dark matter particle mass; comoving (Plummer-equivalent) gravitational softening; maximum physical softening, and a 
brief description. 


Simulation 

L 

(cMpc) 

N 

mb 

(M0) 

^dm 

(Mq) 

Ccom 

(ckpc) 

Cprop 

(pkpc) 

Remarks 

Ref-L100N1504 

100 

2 X 15048 

1.81 

X 

10® 

9.7 X 10® 

2.66 

0.7 

ref. stellar &: ref. AGN feedback 

Ref-L050N0752 

50 

2 X 7523 

1.81 

X 

10® 

9.70 X 10® 

2.66 

0.70 


Ref-L025N0376 

25 

2 X 3763 

1.81 

X 

10® 

9.70 X 10® 

2.66 

0.70 


Ref-L025N0752 

25 

2 X 7523 

2.26 

X 

10® 

1.21 X 10® 

1.33 

0.35 

? ? 

Recal-L025N0752 

25 

2 X 7523 

2.26 

X 

10® 

1.21 X 10® 

1.33 

0.35 

recalibrated stellar & AGN feedback 

NoAGN 

25 

2 X 3763 

1.81 

X 

10® 

9.70 X 10® 

2.66 

0.70 

ref. stellar Sz no AGN feedback 

WeakFB 

25 

2 X 3763 

1.81 

X 

10® 

9.70 X 10® 

2.66 

0.70 

weak stellar & ref. AGN feedback 

StrongFB 

25 

2 X 3763 

1.81 

X 

10® 

9.70 X 10® 

2.66 

0.70 

strong stellar & ref. AGN feedback 

NoFB 

25 

2 X 3763 

1.81 

X 

10® 

9.70 X 10® 

2.66 

0.7 

no feedback using Anarchy SPH 

NoFB-Gadget 

25 

2 X 3763 

1.81 

X 

10® 

9.70 X 10® 

2.66 

0.70 

no feedback using Gadget SPH 

Ref-Gadget 

25 

2 X 3763 

1.81 

X 

10® 

9.70 X 10® 

2.66 

0.70 

Ref-L025N0376 using Gadget SPH 


the average density of the halo equals 200 times the mean 
density of the Universe at any given redshift. The mass con¬ 
tained within that radius is then defined as the halo mass, 
M 200 . The most massive snbstructure in each halo is de¬ 
fined as the central galaxy. The focus of this study is on 
the distribution of gas around bright high-redshift galaxies 
(e.g., LBGs and bright quasars) which are in most cases the 
brightest and most massive objects in their host halos, we 
only consider central galaxies in our analysis. 


2.3 Hi fractions 


For an accurate calculation of the simulated Hi column den¬ 
sities the main ionizing processes that shape the distribution 
of neutral hydrogen must be taken into account. Besides col- 
lisional ionization, which is dominant at high temperatures, 
photoionization by the metagalactic UVB radiation is the 
main contributor to the bulk of hydr ogen ionization on cos - 
mic scales, particularly at z > 1 (e.g.. lRahmati et al.ll2013^ ). 
On smaller scales and close to sources, local radiation could 
be t he dominant source of photoionization (see Appendix m 
and iRahmati et al.l[2013bl ). 

In this work, we use the UVB model of 
iHaardt fc Madaul ll200ll) to account for the mean ionizing 
radiation held from quasars and galaxies. This UVB model 
was also used for calculating radiative heating/cooling 
rates in the hydrodynamical simulations. Moreover, this 
UVB model has been shown to reproduce results consistent 
with the observed Hi column density distribution function 
I Rahmati_et_aljj2013a|) and 2^3 metal absorption lines 
I Aguirre et al. l2008l) . However, we emphasis that both 
observational constraints and model predictions for the 
amplitude (Fuvb) and spectral shape of the photo-ionizing 
background are uncertain by a factor of a few (e.g. . 


Faucher-Giguere et al.l l2008l . l2009l : iHaardt fc Madaul I2OI2I : 


Becker fc Boltor] 20131 ). which could change the Hi column 


density distribution. Where appropriate, we consider the 


impact of those uncertainties on our results by varying the 
UVB model in our calculations. 

At very low Hi column densities, where the gas is highly 
ionized the gas satisfies the so-called optically-thin limit. 
As the Hi column density increases and its corresponding 
optical depth gets close to unity, photon absorptions be¬ 
come more important and eventually the gas can shield it¬ 
self against the incoming flux of ionizing photons. To ac¬ 
count for this self-shielding, we use the same approach as 
we adopted in IRahmati fc Schay 3 ||2014|. Namely, we use 
the htting function presented in IRahmati et akl (l2013al ) for 
calculating the photoionization rate and hence the ioniza¬ 
tion state of hydrogen atoms. This htting function accu¬ 
rately reproduces the result from radiative transfer simula¬ 
tions of the UVB and recombinatio n radiation in cosmolog- 
ical d ensity helds using T RAPHIG (jPawlik fc Schav3l200^ . 
I 2 OIII : iRaicevic et al.ll2013l ) . One can characterize the UVB 
at any given redshift by the hydrogen photoionization rate in 
optically-thin limit, Fuvb.z, and the effective hydrogen ab¬ 
sorption cross-section, (Juui- Then the htting function gives 
the effective photoionization rate, Fphot, as a function of 
density: 


Fphot 

Fuvb.z 


= 0.98 


1 + 


n„ 


UH.SSh 


+0.02 


1 + 




«H,SSh 


( 1 ) 


where is the hydrogen number density and nn.ssh is 
the self-s hielding density threshold predicted by the analytic 
if lSchavd (|2001bh 

2/3 

= 6.73 X 10"®cm“® ' 


model o: 

?1-H,SSh 


2.49 X 10-18 cm-2 


Fuvb.z 
10-12 s-i 


2/3 


( 2 ) 
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Figure 1. CDDF of neutral gas at different redshifts for the 
Ref-L100N1504 EAGLE simulation. The data points represent 
a compilation of various quasar absorption line observations 
at hig h redshifts (i.e., 2: = [ 1.7,5.5]) taken f rom I Peroux et al.l 
ll200a ) wi th 2: = [1.8,3.5], lO^k^ara et al.l ll200'!il) with 2 = 
[1.7,4.5 1, iNoterdaerne et al.l l|2Q09^ with ^ = [2.2,5.5] and 

IProchaska &: Wolfel <l2009h with 2 = [2.2,5.5]. The grey dia¬ 
monds at A^hi > 10^^ cm“^ represent the most recent con- 
straints on the high end of the Hi CDDF which are taken from 
iNoterdaeme et al.l l|2Q12h with (z) = 2.5. The grey star-shaped 
data-p oints at A^hi < lO^^cm”^ are taken from iRudie et ^ 
l|2013l ) with 2 = [2.0,2.8]. The simulation results are in good 
agreement with the observations. 


Combining equations o and © allows us to calcu¬ 
late the equilibrium hydrogen neutral fraction of each SPH 
particle in our hydrodynamical simulations after calculating 
its collisional ionization rate, which depends on the temper¬ 
ature, and its optically-thin recombination rate (i.e.. Case 
Afl which depen d on both the density and temperature (see 
Appendix A2 of fRahmati et aLll2013al l. Since the tempera¬ 
ture of star-forming gas in our simulations is defined by a 
polytropic equation of state that is used to limit the Jeans 
mass, and therefore is not physical, we set the temperature 
of the ISM particles to Tism = 10^ K which is the typical 
temperature of the warm-neutral ISM. 

Noting that the covering fraction of extremely high Hi 
column densities are negligible, we chose not to account 
for the formation of molecular hydrogen wh ich is expected 
to be dominant only at Nm > 10^^ cm~^ l|Schavel f2001cl : 
iKrumholz et al.| l2009l : iRahmati et al.l[2013ah . We also ne- 


^ Note that the impact of recombination radiation is accounted 
for by the fitting function and there is no need to use the Case B 
recombination rate. 
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Figure 2. Cosmic density of Hi as a function of redshift in the 
Ref-L100N1504 EAGLE simulation (solid curve). The shaded 
area around the curve indicates the range covered by all the 
simulations listed in Table ^ (expect models with no feedback). 
The data points represent a compil ation of various q uasar ab¬ 
sorption line observations taken from lRao et al.| <|20 Q6h with 2 = 
[0.11 — 1.65]' [^ochaska et al.l ll2 005h: jProchaska &; Wolfel ll2009l) 
with 2 = [2 .2, 5.5], INoterdaeme et alj ||2012|) with 2 = [2 . 0,3.5 ], 
IZafar et al.l ll 20131) with 2 = [1.5, 5.0] and ICrighton et alj ll2015h . 
The low- redshift comp i lation of data is based on 2 1 cm emission 
stu dies oflZwaan et al.l < |2005l) , iMartin et al.l 1I2OIOI) , iBraunI <|2012|) 
and lPe lhaize et al.l ll2013h at 2 0, stacked 21 cm emission stud¬ 

ies of iRhee et al.l (l2013h at ^ ~ 0.1 — 0.2 and 21 cm intensity 
mapping of lChangeT'^ ll2010l i at z ^ 0.8. 


gleet the impact of radiation from local sources, which is 
thought to become increasingly important for very high H i 
column densities and very close to galaxies (ISchavd l200fil : 
IRahmati et al.|[2013bfi . Noting that the distribution of LLSs 
around the virial radii of galaxies is n ot strongly affected by 
local radiation (see Appendix IbI and IRahmati et al.l l2013bl : 
IShen et al.l l2013l : IRahmati fc SchavelkoilT we postpone a 
treatment of local radiation, which is potentially important 
but requires complex radiative transfer simulations, to fu¬ 
ture work. 

To calculate Hi column densities we use SPH interpo¬ 
lation and project the Hi content of desired regions, which 
can range from the full simulation box to only a small vol¬ 
ume around a galaxy, onto a 2-D grid. We found that using 
a grid with 10, 000^ = 10® pixels for projecting the full box 
of the Ref-L100N1504 simulation results in converged Hi 
covering fractions for Nm < 10^^cm~^, which is the range 
of Hi column densities we study in this work. We use 16 
slices with equal widths for calculating the Hi column den¬ 
sities in the full 100 Mpe simulation box. This enables us 
to calculate Hi column densities as low as Nm ~ 10*^® cm“^ 
without being affected by projection effects. Moreover, this 
choice enables us to calculate the covering fraction of Hi 
around simulated galaxies in analogy with what is done ob- 
servationally, where absorbers are considered to be around 
a galaxy only if their line-of-sight velocity differences from 
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that of the galaxy do not exceed a hxed valu^l- Splitting 
the full 100 cMpc simulation box into 16 slices results in a 
velocity width few times smaller than the typical velocity 
cuts used in observational studies (e.g., AV « 400 km/s at 
z ^ 2 for each slice). Adding together appropriate number of 
slices allows us to efficiently calculate the distribution of Hi 
around galaxies by using velocity cuts compar able to what 
is typically used in observati onal studies (e.g., iRudie et al.l 
1201^ IProchaska et al.ll2013allb[l . 


In previous theoretical studies, the covering fraction 
of Hi was usually measured by considering a finite region 
around galaxies which is o ften normalised to the vir i al ra¬ 
dius of each galaxy ( e .g.. iFa u cher-Giguere fc Kere^ 


Faucher-Gigu^e et al.l 12015 r iFumagalli et al.l l201ll . 


2011 ; 


2014; 


Shen et al. 2013l f . Tdie covering fractions measured using 


this method only account for the gas that is much closer 
to galaxies along the line-of-sight compared to what is 
measured observationally. Observations take into account a 
much longer path-length along the line-of-sight compared 
to the virial radii of galaxies and are therefore not directly 
comparable to the covering fractions reported in previous 
theoretical studies. However, considering a hnite region nor¬ 
malised to the virial radius of each galaxy is a sensible choice 
to study the distribution of gas that has a close physical 
connection to galaxies. For this reason, and also to facili¬ 
tate comparison between our results and previous theoretical 
work, we not only mimic the observed velocity cuts, but also 
measure the covering fraction of Hi systems around galaxies 
by considering only absorbers that are within 2 x r 2 oo from 
each galaxy. 


3 RESULTS 


3.1 Cosmic distribution of Hi 


Before studying the distribution of Hi around simulated 
galaxies, we need to test whether the EAGLE simulation 
reproduces the observed cosmic distribution of strong Hi 
systems. If the simulation fails to satisfy the existing con¬ 
straints on the cosmic abundance of Hi absorbers, then it 
is hard to trust the predictions drawn from it about the 
distribution of Hi close to galaxies. While performing this 
consistency check is not straightforward for studies that use 
zoom simulations, using a cosmological simulation enables 
us to do so. 

The cosmic Hi distribution is often quantified as the 
Hi column density distribution function (CDDF). The Hi 
GDDE, f{Nm,z), is conventionally defined as the num¬ 
ber of absorbers per unit column density, dNm, per unit 
absorption length, dX = dz {Ho/H{z)){l + z)'^, and is 
measured observationally by searching for Hi absorbers 


2002 

. I 2 OI 3 I: iPeroux et al.l I 2 OO 5 I: lO’Meara et al. 

2007, 

2013 

: Noterdaeme et al. 20091. 20121: IProchaska et al. 

2009; 

Prochaska & Wolfe l2009l: Zafar et al.l 20131: Rudie et al. 


^ The typical velocity differe nce cut often used in observational 
studie s is > ±1000 km/s (e.g.. [Rudie et aI1l2012l : IProchaska et all 
l2013allbh . 


the Ref-Ll DON1504 simulation with a compilation of ob¬ 
servational result^. The predicted Hi GDDE is shown using 
different colours and line styles for different redshifts ranging 
from 2 ; = 1 to 5. For comparison, the grey data points show 
the observed Hi CDDF for Nm < 10^^ cm“^ at 2 ~ 2 — 3 
from iRudie et al.l (l2013f ). and at higher Ahi a compila¬ 
tion containing various observations spann i ng the redshift 
range of 1.7 < 2 < 5.5 ( Perouxet_aL 200^ ^MearaetaJj 


1200^ : iNoterdaeme et ah l2009l : IProchaska fc Wolfel l2009l ~l. 


The most recent measurements of the Hi CDDF at very 
high A^hi and an aver age redshift of ( 2 ) = 2.5 from 
INoterdaeme et ^ ll2012fl are also shown using dark grey di¬ 
amonds. 

Fig. [ 1 ] shows that there is good agreement between the 
predicted Hi CDDFs and observations for strong Hi ab¬ 
sorbers (Nm > 10^^cm~^), similar to what wa s found for 
OWLS llAltav et al.ll20lI] : lRahmati et al.ll2013al l. The weak 
evolution of t he high end of the Hi C DDF that we reported 
for OWLS in iRahmati et al.l ll2013ah is also evident . How - 
ever, we note that the measurements of IRudie et al.l (l2013fl 
for the CDDF are slightly under-produced by our simulation 
which indicates the need to use a lower hydrogen photoion¬ 
ization rate at 2 = 2.5 compared to what our fiducial UVB 
model implies (see Appendix o. It should also be noted 
that because we do not correct the simulation for H 2 , the 
agreement at Nm ^ 10^^ cm“^ may be fortuitous. 

The cosmic Hi density, Dhi, which is defined as the 
mean Hi density divided by the critical density, pcrit, can be 
calculated either by estimating the observed Hi mass density 
through Hi 21-cm emission at low redshifts or by integrating 
the Hi CDDF of absorbers at high redshifts: 


flni = 


Homn 

Cpcrit 


f 


Nmf(Nm, 2)dAHi, 


(3) 


where Ho — 100 h km s“^ Mpc“^ is the Hubble constant, 
mn is the mass of a hydrogen atom, c is the speed of light 
and Pcrit = 1.89 x 10~^^h^g cm“®. 

The predicted evolution of the cosmic Hi density for 
the Ref-L100N1504 EAGLE simulation is shown with the 
solid curve in Fig. [51 The shaded area around the curve 
shows the range covered by the other feedback enabled sim¬ 
ulations we use in this work (see Table [T|. For compari¬ 
son, a compilation of various observational measurements 
are over plotted with different symbols. The high -2 (2 > 1) 
measurements of the flni are often based on the observed 
abun d ance of DLAfl(e.g., iRao et ahllio o^: IProchaska et aP 


2 OO 5 I: Prochask a ifeWolfd I 2009 I : Noterdaeme et akT 20121 : 
Zafar et al ■II2013I) . The low-redshift measurements are gener¬ 
ally based on measuring the Hi mass using 21-cm emission 
and often involve adopting non-trivial assumptions about 
the Hi gas fraction of the full galaxy population to derive 
the c o smic Hi density (e.g ..lZwa an et al.ll2005l : iMartin et al 


20K; IChang et al.l I 2 OI 0 I : IPelhaize et al.l 1201^ 


Rhee et al 


2013ll . The comparison between the predicted and observed 
Dhi shows good agreement, particularly at the redshifts of 


® We apply appropriate corrections for the different cosmological 
parameters adopted in different studies, converting all results to 
the Planck cosmology used in EAGLE. 

^ Note that due to the shape of the Hi CDDF, the cosmic Hi 
density is very sensitive to the abundance of DLAs and that the 
contribution of lower Nm systems is negligible. 
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Logio [Nh, (cm"^)] 


Figure 3. The simulated Hi column density distribution around randomly selected massive galaxies at z = 3. Top, middle and bottom 
rows show galaxies with M 200 ~ 10^^, ~ and 10^^ Mq, respectively. The columns in each row show a single galaxy as seen from 

three different orthogonal angles. Blue circles are centred on galaxies and show the virial radii (i.e., r20o)- Each panel shows a 1 x 1 pMpc^ 
region with the same projected depth. The covering fraction of LLSs (i.e., A^hi > cm“^) with impact parameters less than r2oo is 

indicated in the top-right of each panel. At z = 3, LLSs (green regions) form filamentary structures and their distribution varies strongly 
from galaxy to galaxy, and with the viewing angle. The typical covering fraction of LLSs within r2oo does not vary strongly with halo 
mass at a given redshift. 
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Figure 4. The same as Figure [3] but for galaxies with M 200 ~ Mq at redshifts 2 : = 4, 2 : = 3 and z = 2 (from top to bottom, 

respectively). The columns in each row show a single galaxy as seen from three different orthogonal angles. LLSs (green regions) form 
filamentary structures at all redshifts, which become less prominent by decreasing redshift. The typical covering fraction of LLSs within 
^200 evolves rapidly with redshift. 
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interest here, z > 1. While the cosmic Hi density remains 
nearly constant from 2~6toj2~2 — 3, it declines to¬ 
wards lower redshifts. This decline, which is also evident in 
observed evolution of the flni, is not reproduced in some pre¬ 
vious theoretical studies le.g.. lDave et al.ll2013l : iLaeos et alJ 
|20141 . The cosmic Hi density in the Ref-L100N1504 simu¬ 
lation seems to drop faster than observed while the range of 
predictions obtained by varying the resolution and/or box- 
size of the simulation (shaded region around the solid curve 
in Fig. [2J is fully consistent with the observational measure¬ 
ments at 2 « 0. Moreover, as we mentioned above, it is 
important to note that measurements of flni at low redshift 
often involve strong assumptions about the Hi mass frac¬ 
tions of all galaxies and are therefore not as robust as the 
direct measurements at higher redshifts. 

Having shown that the observed cosmic distribution of 
Hi is reproduced reasonably well, we can use the simulation 
with more confidence to study the Hi distribution around 
galaxies. 

3.2 Covering fraction of Lyman Limit Systems 
inside halos 

Examples of the distribution of Hi around simulated galax¬ 
ies at 2 = 3 are shown in Figs. [3] In this figure, the 
coloured maps show the Hi column density distributions 
in 1 X 1 pMpc^ regions centred on galaxies with M 200 = 
10^^ — 10^^ Mq. The virial radius, r 2 oo, of each galaxy is 
shown with a blue circle and each galaxy map is shown 
using three different orthogonal projections. A significant 
fraction of the area within the virial radii of massive galax¬ 
ies is covered with LLSs (i.e., A^hi > cm“^) which 

have highly inhomogeneous distributions with often form 
filamentary structures. As a result, the fraction of the area 
covered by LLSs inside the virial radius, which is indicated 
in the top-right corner of each panel, varies depending on 
the point of view and from one galaxy to another. However, 
the average LLSs covering fraction does not depend strongly 
on the halo mass. On the other hand, as|3]shows, the typical 
LLSs covering fraction decreases significantly from 2 = 4 to 
2 . 

To quantify the distribution of Hi around galaxies, the 
covering fraction of LLSs within the virial radius, f<r 2 oo, 
is defined as the probability of finding systems with A^hi > 
cm“^ with impact parameters smaller than the virial 
radius, r 2 oo and with line-of-sight distances from the galaxy 
shorter than a specific value comparable to the virial radius. 
Equivalently, /<r 2 oo defined as the fraction of the 

surface area within R < r 2 oo that is covered by LLSs after 
projecting the gas distribution within a specific line-of-sight 
distance from the galaxy onto a 2-D plane. We calculate 
this quantity for each galaxy by projecting the Hi within a 
slice with 4 x r 2 oo thickness centred on the galaxies redshiflQ 
and repeating the same calculation for projections along 3 
different orthogonal directions. Then we calculate /<r 2 oo by 
measuring the fraction of the surface contained within the 
r 2 oo that is covered by LLSs. 

® We us ed this definition to be consistent with previous stud¬ 
ies (e.g.. I^magalli et al.|[?014^ . We note, however, that choosing 
thinner or thicker slices does not change our results noticeably. 


The predicted f<r 2 oo different redshifts is shown in 
Fig. [5]as a function of halo mass, M200, in the left panel, and 
as a function of specific star formation rate, M*/M*, in the 
right panel. Each solid curve and the shaded area around 
it show, respectively, the median covering fraction and the 
corresponding 15-85 percentiles for the Ref-L100N1504 sim¬ 
ulation. In this figure, red, orange, green and blue curves 
show 2 = 4, 3, 2 and 1, respectively. For massive haloes 
with M200 ^ 10^^ Mq, f<r 2 oo do6s not depend strongly on 
halo mass. For less massive haloes (i.e., M200 10^^ Mq) on 

the other hand, the covering fraction of LLSs increases more 
strongly (though still relatively weakly) with halo mass and 
the slope of this relation increases with redshift. 

While the dependence of /<r 2 oo halo mass is rather 
weak, the covering fractions increase strongly with redshift. 
This result is consistent with halos containing increasingly 
higher gas fractions at higher redshifts as a result of increas¬ 
ing rates of cold accretion and the higher mean density of 
the Universe. As we will show in this weak halo mass 
dependence enables us to characterise the distribution of ab¬ 
sorbers over a wide range of halo mass as a function of red¬ 
shift and radius relative to the virial radius. The strong cor¬ 
relation between the covering fraction of LLSs and redshift 
has important consequences for the interpretation of obser¬ 
vations because the observed sample often contain galaxies 
with a wide range of redshifts. For instance, the average 
probability of finding LLSs within a given distance from 
galaxies does not represent the covering fraction of LLSs 
at the typical (e.g., mean) redshift of the galaxy sample, be¬ 
cause higher redshift galaxies have larger contributions to 
the average covering fraction of the population. As we will 
discuss in m together with other biases like the wide range 
of halo masses represented by observational samples, this 
issue can explain the large cover ing fractions derived fro m 
some observational samples le.g.. IfTochaska et al.ll^l3ljj . 

The right panel of Fig.[5]shows that /<r 2 oo does not de¬ 
pend strongly on the specific star formation rate of galaxie^. 
Only galaxies with M200 > Mq are shown in the right 

panel of Fig. O but our experiments show that a narrow 
range of halo masses only strengthens the independence of 
f<r 2 oo from the specific star formation rate. This trend thus 
suggests that the covering fraction of LLSs does not depend 
strongly on the transient variations in the star formation ac¬ 
tivity of galaxies, and is set by their average star formation 
activity and the large-scale distribution of gas around them. 

As the shaded areas around the median curves in Fig. 
[S] illustrate, there is significant scatter in the predicted cov¬ 
ering fraction at any given mass and specific star formation 
rate. This scatter is larger at higher redshifts where the typ¬ 
ical covering fractions are also larger. We note that the cov¬ 
ering fraction of a single simulated galaxy can change from 
one projection axis to another by a factor close to the typical 
scatter fo r its mass range (see F ig. [3]l, which is consistent 
with what lFumaealli et al.l ||2014|') found. This, together with 
the lack of strong dependence between the covering frac¬ 
tion and specific star formation rate, suggests that most of 
the scatter shown in Fig. [5] can be attributed to the highly 


® Note that we neglect the impact of local sources on /<r 2 oo- 
local sources were to change the f<r 2 oo significantly, then they 
could introduce a dependence on the specific star formation rate. 
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Figure 5. Cumulative covering fraction of LLSs within r2oo a function of halo mass, M 200 (left) and specific star formation rate, 
M*/M* (right). Curves from top to bottom show redshifts z = 4 (red), 3 (orange), 2 (green) and 1 (blue). Solid curves show the median 
covering fractions for the Ref-L100N1504 simulation and the shaded areas around them indicate the 15 — 85 percentiles. Galaxies are 
shown with individual data points instead of curves in bins that contain fewer than 10 galaxies. About 35,000 galaxies (each in 3 different 
orientations) were used to make this figure. The covering fraction increases strongly with redshift for all halo masses. It also increases 
with halo mass but this dependence is weak for massive objects. There is no strong correlation between the covering fraction of LLSs 
and the specific star formation rate. 


inhomogeneous and filamentary distribution of Hi around 
galaxies. 

Although f<r 2 oo widely used to quantify the distri¬ 
bution of Hi around galaxies, both in theoretical and ob¬ 
servational studies, one should note that the virial radius is 
not a directly observable quantity. Moreover, as mentioned 
above, the virial radius of a sample of galaxies with a wide 
range of different characteristics (e.g., mass, redshift) is not 
well defined. For this reason we opted not to compare the 
covering fractions shown in Fig. [5] with tho s e reported in ob¬ 
servat ional studies (e.g.. lRudie et al']l2012l : [Prochaska et al.l 
l2013all^ . Instead, we compare to observations after match¬ 
ing the redshift and halo mass distribution of observed and 
simulated samples in Sec [I] 

3.3 Covering fraction profiles 

In the previous section we studied the cumulative covering 
fraction of LLSs within the virial radius (i.e., f<r 2 oo)- How¬ 
ever, more information is embedded in the differential cov¬ 
ering fraction profile of Hi absorbers with different column 
densities as a function of impact parameter from galaxies. 
We define the differential covering fraction in a given impact 
parameter bin as, 

Aaba 

fcav{R) = fcov{Ri < R < Ri+l) = „ 2\ ’ 

where Aabs is the area covered by absorbers (e.g., LLSs) 
within a radial bin defined by the impact parameters Ri and 
Ri+i, assuming that Ri+i > Ri > 0 . Note that throughout 
this work we reserve fcov{R) the differential covering frac¬ 
tion and /cov(< R) for the cumulative covering fraction (e.g.. 


Figs. 151 and llOD . For instance, we denote the cumulative cov¬ 
ering fraction within an impact parameter equal to the virial 
radius by /<r 2 oo = /cov(0 < R < r2oo) as shown in Fig. [S] 

The top-left panel in Fig. [6] shows the predicted mean 
differential covering fraction of LLSs as a function of im¬ 
pact parameter for five different halo mass bins in the Ref- 
LIOONI5O4 simulat ion at z = 2.5. In analo g with observa¬ 
tional studies (e.g.. llTochaska et al.ll2013alra) . around each 
galaxy a velocity window of AV = 3150km/s (i.e., the al¬ 
lowed velocity difference between absorbers and galaxies is 
< ±1575 km/s) is adopted for calculating its covering frac¬ 
tion, but we note that increasing or decreasing the allowed 
velocity width by a factor of a few does not change the re¬ 
sults for R < r 2 oo (see Fig. lCljl . The five mass bins shown in 
the top-left panel of Fig. [6] have similar LLSs covering frac¬ 
tions at the outermost impact parameters, but they vary 
strongly with halo mass close to galaxies. As the middle left 
and bottom left panels in Fig. [6] show, the same qualitative 
trend holds for sub DLAs (i.e., Ahi > 10^^ cm“^) and DLAs. 
However, by increasing the Hi column density of absorbers, 
the covering fraction at fixed impact parameters decreases. 
Despite the sensitivity of the covering fractions to the halo 
mass, it seems that the shapes of the curves are very similar 
for halo masses M 200 ^ 10^^ Mq which suggests that they 
can be matched by a re-scaling to account for differences in 
the virial radii of the halos. 

To show that the covering fraction profiles are nearly 
scale invariant, we normalize the impact parameters to the 
virial radii of the halos in the panels of the right side of 
Fig. The good agreement between the three curves that 
represent M 200 > 10*^^ M© halos shows that the covering 
fraction profiles of LLSs, sub DLAs and DLAs around those 
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Figure 6. Profiles showing the mean differential covering fraction of LLSs (top), sub DLAs (middle) and DLAs (bottom) as a function of 
impact parameter (left) and normalised impact parameter (right) for different halo mass bins in the Ref-L100N1504 EAGLE simulation 
at 2 = 2.5 using AV = 3150 km/s. In the left panel, the virial radius corresponding to each halo mass bin is indicated with a cross. As the 
left panels show, all covering fractions depend strongly on halo mass at fixed impact parameters, particularly very close to galaxies. The 
right panels show, however, that only a weak mass dependence remains in the covering fraction profiles of halos with M 200 ^ 10 ^^ Mq, 
after normalizing the impact parameters to the virial radii. This suggests that the shape of the covering fraction profiles of LLSs, sub 
DLAs and DLAs is similar in all halos with M 200 ^ 10^^ Mq and the mass dependence of the covering fractions at fixed physical impact 
parameters stems mainly from the differences in the halo sizes. 
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Table 2. The best-fit values for the free parameters of the fitting function, equations [5]&:|6l for the predicted normalised covering fraction 
of LLSs, sub DLAs and DLAs for halos with M 200 > 10^^ Mq at redshifts 2 < 4. Note that parameter C is sensitive to the chosen 
velocity width and is reported for two velocity width AV = 3000 and 1500 km/s. The performance of the fitting function (thin solid 
curves) is shown in Figure[7]for AV = 3000km/s. 


Absorbers 

A 

B 

C 

C 

a. 




AH = 3000 km/s 

AH = 1500 km/s 


LLS (Am > 10i'^-2 cm2) 

0.100 

1.90 

0.21 

0.15 

2.1 

sDLA (Am > lOi'^ Ocm^) 

0.060 

1.83 

0.07 

0.05 

2.0 

DLA (Am > 1020'3cm2) 

0.035 

1.73 

0.02 

0.015 

1.8 


halos are self-similar with a characteristic scale length very 
close to the virial radius. This is the main reason behind the 
very weak dependence of the total LLS covering fraction 
within r 2 oo (/<r 2 oo) halo mass for galaxies with M 200 ^ 
10^^ Mq (see Fig.[S]|. As the curves for the two lowest mass 
bins in the right panels of Fig. |6]show, the scale-invariance 
of the covering fraction profiles breaks down for galaxies 
with M 200 < 10^^ Mq, where the total covering fraction of 
absorbers within r 2 oo also slowly decreases with decreasing 
halo mass (see Fig. O. 

The scale-invariance of the distribution of strong Hi ab¬ 
sorbers around massive halos allows us to calculate the typ¬ 
ical normalised covering fraction profiles that characterise 
the distribution of LLSs, sub DLAs and DLAs around ha¬ 
los with M 200 g/ 10*^^ Mq at any given redshiflQ. Based on 
the strong evolution of the total Hi covering fraction in¬ 
side halos (see Fig. 0), it is expected that the normalised 
covering fraction profiles also evolve strongly with redshift. 
To illustrate this, we show the normalised differential cover¬ 
ing fraction for LLSs, sub DLAs and DLAs in, respectively, 
the top, middle and bottom panels of Fig. [7] The differ¬ 
ent curves in each panel indicate different redshifts where 
long-dashed (red), dashed (orange), dot-dashed (green) and 
dotted curves show 2 = 4, 3, 2 & 1, respectively. Note that 
a line-of-sight velocity window of width AH = 3000 km/s 
is adopted for calculating the covering fractions shown in 
this figure. To illustrate the typical intrinsic scatter in the 
covering fraction profiles, the 15-85 percentiles of the cov¬ 
ering fractions at 2 = 4 are indicated by the shaded areas 
around the long-dashed red curves. The normalised cover¬ 
ing fraction profiles of all strong Hi absorbers indeed evolve 
strongly. 

Defining x = R/r 200 as the normalised impact param¬ 
eter, the covering fraction of Hi absorbers around galaxies 
with a given virial radius, r 2 oo, at redshift 2 can be fit with 
the following function: 

10^, (5) 

where is a redshift-dependent characteristic length scale 
that is given by: 

L,=A (6) 

and A, B, C and a are free parameters that vary for LLSs, 


^ The Ref-L100N1504 EAGLE simulation contains 95, 345, 824 
and 1436 halos with M 200 ^ 10^^ Mq at redshifts 2 = 4, 3, 2 and 
1, respectively. 


sub DLAs and DLAs. Based on this fitting function, the cov¬ 
ering fraction approaches unity very close to galaxies where 
X Lz and at very high redshifts where 2 —>■ 00. However, 
note that the latter is the case only if B < 10*^^®. Far from 
galaxies where x ^ Lz, on the other hand, the covering 

z -4 

fraction approaches the asymptotic value of C 10 3 . 

Table [ 2 ] lists the best-fit values for the free parameters. 
As shown in Appendix El the covering fraction of absorbers 
at relatively large impact parameters {R > r 2 oo) is sensitive 
to the size of the line-of-sight velocity window that is used 
for associating Hi absorbers and galaxies. As a result, pa¬ 
rameter C in the fitting function of equation m is sensitive 
to the adopted velocity width. Therefore, we report two sets 
of best-fit values of C where each set corresponds to either 
AH = 3000 km/s or AH = 1500 km/s. 

As the thin solid curves in Fig. [7] show, the predicted 
normalised covering fraction profiles are all closely matched 
by the values obtained from equation © for appropriate 
choices of the free parameters that are listed in Table (2] We 
note that the differences between the fitting formula and the 
simulation are much smaller than the typical intrinsic scatter 
in the covering fractions, which are shown in Fig.Qfor 2 = 4 
by the shaded areas around the long-dashed red curves. Note 
that variations in the assumed UVB radiation can change 
the covering fraction of LLSs. For instance, reducing the 
UVB photoionization rate by a factor of 3 results in LLS 
covering fractions that are higher by 0.1. However, such 
moderate changes in the UVB do not significantly affect the 
distribution of highly self-shielded stronger absorbers such 
as sub DLAs and DLAfl. 

The values of the free parameters that control the em¬ 
pirical fitting function that we introduced above are phys¬ 
ically meaningful. For instance, Lz could be regarded as a 
typical projected distance between absorbers and galaxies. 
Taking the values of A and B from Tabled the implied typi¬ 
cal projected distances between LLSs, sub DLAs and DLAs 
and their host galaxies at 2 « 3 are respectively « r 2 oo, 
« 0.5r2oo and » 0.2r2oo, which is i n exc ellent agreement 
the predictions of iR.ahmati fc Schavd (l2014l l from the OWLS 
simulations (see the right panel of their Fig. 3). Moreover, 
equation implies that the typical normalized impact pa¬ 
rameter of strong Hi absorbers increases exponentially with 
redshift. The best-fit values for B, which controls the rate 
of this change, imply that the impact parameters of LLSs 
evolve slightly faster than those of DLAs. Moreover, the 

® It is important to keep in mind that very close to galax¬ 
ies (R ^ r20o) the predicted covering fractions may be over¬ 
estimates because we neglect local sources of ionizing radiation. 
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Figure 7. Normalised differential covering fraction profiles of 
Hi absorbers around galaxies with M 200 ^ 10^^ Mq in the 
Ref-L100N1504 EAGLE simulation at different redshifts using 
AV = 3000 km/s. The top, middle and bottom panels show the 
covering fractions of LLSs, sub DLAs and DLAs, respectively. 
The red (long-dashed), orange (dashed), green (dot-dashed) and 
blue (dotted) curves show the results at z = 4, 3, 2 and 1, re¬ 
spectively, while the thin solid curves show the fitting function 
given by equation [5] and [6] The covering fractions for stronger Hi 
absorbers are lower and their profiles have shorter characteristic 
scale lengths. The covering fractions of all absorbers drop rapidly 
with decreasing redshift at all impact parameters. 


higher a value for LLSs indicates that their covering frac¬ 
tions drop faster than those of DLAs with increasing nor¬ 
malised impact parameters. 

The right-most term in equation m is related to the 
contribution of the background absorbers outside of the 
halo. Given the steep Hi CDDF (see Fig.[TJ, there are many 
more LLSs than DLAs. Given the Hxed line-of-sight velocity 
cut imposed to obtain the covering fraction of absorbers, one 
would expect the ratio between the covering fraction of ab¬ 
sorbers at large virial radii to follow the ratio between their 
cosmic abundances. As the best-fit values for the C param¬ 
eter imply, this is indeed the case (e.g., LLSs are ~ 10 times 
more frequent than DLAs at all redshifts). 

We note that local sources of ionising radiation, which 
we have ignored, may cause the fitting function to under¬ 
predict the covering fraction of strong Hi absorbers close 
to galaxies (see Appendix [Q. We emphasise that we only 
considered the redshift range 1 < 2 < 4 when deriving our 
htting function. While the same function produces a reason¬ 
able match to the simulated normalised Hi covering fractions 
at z = 5 for R < r 2 oo, it predicts LLS covering fractions 
that are « 10% too high for larger impact parameters. This 
difference is smaller than (but comparable to) the intrinsic 
scatter (i.e., 15-85 percentiles) in the simulated normalized 
profiles. 


4 COMPARISON WITH OBSERVATIONS 

Recent observations found large coverin g fractions o f LLSs 
close to massive halo s at z ^ 2 IjRudie et al.l l2012l : 

IProchaska et al.l I2ni3allbh which implies the existence of 
a large reservoir of neutral and hence relatively cold gas 
around massive a ~ 2 galaxies. However ^ recent simula- 
tions of ~ 10^^'^ MfTi halos at g ~ 2 by iFumagalli et al.l 
ll2014h and iFaucher-Giguere et al.l ll2015h resulted in LLSs 
co vering fractions t hat ar e much smaller than observed 
bv IProchaska et al.l (I2ni3bh . This discrepancy may indicate 
that our current theoretical understanding of galaxy forma¬ 
tion and evolution is inadequate. 

The EAGLE simulation is ideal to revisit this problem 
since its volume is sufficiently large to contain a large num¬ 
ber of halos with M 200 ^ 10^^ M© at a ~ 2 — 3 without 
compromising the resolution needed to reproduce the ob¬ 
served cosmic distribution of Hi (see 1)3.111 . Because of the 
large intrinsic scatter in the covering fractions, a large sam¬ 
ple of simulated galaxies is required to constrain the average 
covering fraction at any given mass. Moreover, thanks to ef- 
hcient stellar and AGN feedback, EAGLE reproduces the 
basic observed characteris tics of galaxies ove r wide ranges 
of mass and redshift ('S15: llffirlong et aLll2014h . In addition, 
using a large cosmological simulation enables us to calculate 
the Hi covering fractions for a line-of-sight path length com¬ 
parable to what is typically used in observational studies. 

In the following, we compare EAGLEs predictions with 
the observational constraints on the distribution of H i 
around mass ive star-forming galaxi es llRudie et al.l l2012fl 
and quasars (IProchaska et al]l2013bh at 2 : ~ 2. 

4.1 Hi distribution around quasars at 2 ~ 2 

IProchaska et al.l (I2013al f9l observed the distribution of Hi 
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Figure 8. Cumulative redshift distribution of the for eground 
quasars from the QPQ6 sample llProchaska et aill2013bh . Differ¬ 
ent curves show the redshift distribution of foreground quasars in 
different impact parameter bins, identical to those used in Fig. [9] 
Most quasars are at large impact parameters and have a wide 
range of redshifts, extending up to 2 : > 4. More than 30% of 
quasars with R > 100 pkpc have redshifts 2 > 2.5. Given the 
rapid increase of the Hi covering fraction with redshift that we 
found in this work, those quasars can strongly bias the estimated 
covering fraction of LLSs, particularly at large impact parameters. 
Moreover, for a magnitude-limited survey, higher-redshift quasars 
are likely to reside in more massive halos which may further bias 
the estimated covering fractions at fixed impact parameters. 


around bright quasars at 2 : 2. Clustering analysis indi¬ 

cate that those quasars which t ypically reside in massive 
halos with M 200 M© ( White et al.l . Using 

background quasars to probe the distribution of gas in ab¬ 
sorption around a large s ample of foreground qua sars (here¬ 
after the QPQ6 samnleh lProchaska et all (l2013lJ l measured 
the covering fraction of LLSs in different impact parameter 
bins. Adopting a fixed typical halo mass of 10^^ ® Mq at 
z ^ 2, and consequently a fixed vi rial radius of 160 kpc for 
all quasars in the QPQ6 sample, IProchaska et al.l (l2013al l 
concluded that more than > 60% of the area within the 
virial radii of halos with masses around 10^^ ® Mq at a ~ 2 
is cov ered by LLSs. For thei r covering fraction measure¬ 
ment, IProchaska et al.l (1201331 1 adopted a velocity width of 
AV = 3000 km/s to associate absorbers to quasars. 

For a proper c omparison between our simulation and 
the observations of IProchaska et ^ (I2013al f9l. it is impor¬ 
tant to note that the foreground quasars in the QPQ6 sam¬ 
ple have a relatively wide range of bolometric luminosities. 
Since quasars are variable, a one-to-one relation between 
quasar luminosity and halo mass is not expected. However, 
as a rough estimate, if we assume a positive correlation be¬ 
tween the halo mass of galaxies and the bolometric lumi¬ 
nosity of the quasars, we could conclude that quasars in 
the QPQ6 sample represent a relatively wide range of halo 
masses some of which could well exceed the typical halo 
mass of M200 10^^'® Mq. If the host halos of the observed 
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Figure 9. Predicted and observed differential Hi covering frac¬ 
tions around quasars at z 'Ri 2. The data points with error-bars 
show the observations of IProchaska et al.l ll2013bl) for a sample of 
quasar at {z) = 2.3. Predicted mean covering fractions for halos 
with M 200 ^ 1 O^^'^M 0 in the Ref-L100N1504 EAGLE simu¬ 
lation are shown with dark-colored regions which indicate the 
systematic uncertainty in the mean due to uncertainties in the 
hydrogen photoionization rate of the UVB and the redshift range 
of the observed quasars (see the text). The light shaded areas 
indicate the 15 — 85 percentiles for the scatter due to object-to- 
object variation. The top, middle and bottom panels show the 
results for LLSs, sub DLAs and DLAs, r espectively. The s quare s 
in the top panel show predictions from l^magalli et al.l || 203 )- 
The observed covering fractions agree with the EAGLE predic¬ 
tions. 
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quasars indeed have a range of masses, then the adopted 
fixed virial radius of 160 pkpc for calculating the total cov¬ 
ering fraction of LLSs inside the virial radii of halos with 


M 200 


10 " 


Mq at 2 ; ~ 2 could result in an overesti¬ 


mate of the true covering fraction due to the non-negligible 
contribution of halos with M 200 > 10^^ ® Mq in the QPQ6 
sample. 

It is also important to note that the foreground quasars 
in the QPQ6 sample are not all at the same redshift. In ad¬ 
dition, as Fig. [8] shows, there is a systematic trend between 
the typical redshift of foreground quasars and the impact 
parameter at which their gas content is measured. As a re¬ 
sult, the typical redshift of quasars for the smallest impact 
parameters 10 — 100 pkpc) is 2 ~ 2 but it increases to 
2 « 2.5 for impact parameters ~ 1 pMpc. The high -2 tail of 
the distribution for quasars with impact parameters R > 100 
pkpc is quite extended and more than 30% of them have 
2 > 2.5. Given our finding that the Hi distribution around 
galaxies at 2 2 — 3 evolves rapidly, this systematic bias 

should be taken into account when comparing simulations 
with observations. 

In addition, the Hi distribution is sensitive to the in¬ 
tensity of the UVB radiation. However, observational con¬ 
straints on the intensity of the UVB at 2 < 2 < 6 
are model dependent an d uncertain by a factor of a few 
(e.g., iBaitlik et al.l 198^: Rauch et al. 1997: Bolton et ^ 


200^: Faucher- Giguere et al. 200^ : Calverlev et al.l 201ll : 

Becker fc Boltonll2013h . UVB models are also uncertain due 


to the various assumptions they need to adopt (e.g., the 
escape fraction of ionizing photons into the intergalactic 
medium, mean-free-path of ionizing photons, abundance of 
faint sources) and differ from e ach other by a factor of a few 
(e.g., Haj3rdt_fcMjjmJ I 2001 I : iFaucher-Giguere et al.l I 2 OO 9 I : 
iHaardt fc MadaiJ 2012 ). While the inten sity of our fiducial 
UVB model fi.e.. iHaardt fc Madaul[200ll l is well within the 
range of the most recent estimates (e.g., iBecker fc Bolto^ 
l2013l'l . intensities lower by up to a factor of ~ 3 are con- 
sistent with some observations/models at 2 ^ 2 — 3 (e.g ., 
iFaucher-Giguere et akllioo^ . I 2009 I : IHaardt fc Madaull201^ . 
and would further improve the agreement between EA¬ 
GLE and the observed Hi column density distribution below 
Vhi « lO^^cm”^ (see Fig. [T]). 

It is necessary to take the aforementioned considera¬ 
tions into account when comparing simulations and obser¬ 
vations. To do this, we calculate the covering fraction of 
LLSs around simulated galaxies with M 200 > 10^^'® M© in 
the Ref-L100N1504 simulation at 2 = 2.2 and 3, resulting 
in 116 and 39 halojj, respectively, with a median mass of 
M 200 = 10^^ ® MqEj The two selected redshifts bracket the 
range of redshifts that is represented by the QPQ6 sam- 


® Although the number of simulated halos we use is less than the 
observed number (155 simulated halos vs. 646 observed quasars), 
we use ~ 10® sight-lines per simulated object to calculate the 
covering fraction profiles. In other words, we use ~ 10^ sight-lines 
to calculate the predicted Hi distributions that ar e shown in Fig.[^ 
while o nly 600 observed sight-lines are used in iProchaska et aid 
ll2013bl) . 

We note that due to steepness of the mass function around the 
halo mass ranges of interest for our analysis, most selected halos 
have masses close to the lower mass limit we imposed in selecting 
them. Using a similar argument, a small fraction of the observed 


pie. As shown in Appendix El varying the UVB model 
changes the resulting Hi covering fraction. It is therefore 
important to include also the uncertainties in the amplitude 
of the UVB photoionization rate. To do this, we calculate 
the Hi distributions usin g both o ur fiducial UVB model o f 
IHaardt fc Madaiil (1200 il l and the IHaardt fc MadarJ (l2012f l 
UVB model. Noting that at 2 « 3 the latter yields a hy¬ 
drogen photoionization rate ~ 3 times weaker than for our 
fiducial UVB model, we co nsider the 2 = 3 Hi coveri ng frac¬ 
tions calculated using the IHaardt fc Madaul (I 2 OI 2 I I model 
upper limits on the predictions. Given the steep evolution 
in the Hi covering fractions from 2 = 3 to 2, we use the simu¬ 
lation results at 2 = 2 that use our fiducial FHaardt fc Madaul 
1 I 2 OOIII UVB model as lower limits for the predictions. Then, 
we calculate the covering fraction of LLSs, sub DLAs and 
DLAs for each h alo using impact param eter bins identical 
to the analysis of IProchaska et al] (l2013bl l . We use a line-of- 
sight velocity window of AV = 3000 and 3400 km/s around 
each galaxy at 2 = 2.2 and 2 = 3 for calculating the covering 
fractions to mimic closely what is done observationall^F^. 

The predicted covering fractions of Hi absorbers are 
shown in Fig. [9] for LLSs (top panel), sub DLAs (middle 
panel) and DLAs (bottom panel). The upper and lower edges 
of the dark-colored areas in each panel show the lower and 
upper limits of our predicted mean covering fractions ob¬ 
tained by applying the fiducial UVB model to 2 = 2.2 halos 
and a 3 times weaker UVB model to 2 = 3 halos, respec¬ 
tively. The shaded areas around the dark regions, which are 
shown using light colors, indicate the regions enclosed be¬ 
tween the 15 percentiles of the lower limit for covering frac¬ 
tion (i.e., at 2 = 2.2 and using the fiducial UVB model) and 
the 85 percentile of the upper limit for the covering fraction 
(i.e., at 2 = 3 and using the weaker UVB model) in each im¬ 
pact parameter bin. In other words, the dark regions show 
how much variation is expected in the predicted covering 
fractions due to systematic effects caused by the redshift 
distribution of the quasar sample and the photoionization 
rate of the UVB, and the light-colored areas around the dark 
regions show the predicted la scatter (15 — 85 percentiles) 
around the mean due to object-to-object variations in the 
covering fraction within the sample of simulated halos. 

Grey diamond s with error-bars show the observations 
of IProchaska et ^ ll2013bl l where the horizontal error bars 
show the impact parameter bins and the vertical error bars 
show only the Icr statistical uncertainty. Comparing the ob¬ 
served data points with the predicted results shows overall 
good agreement for absorbers with different Hi column den¬ 
sities. The agreement is particularly good for larger impact 
parameters (> 60 pkpc) despite the fact that the obser¬ 
vational error bars are smallest there owing to the larger 
number of quasar pairs. 

For impact parameters in the range 30-60 pkpc we ap¬ 
pear to predict too high covering fractions and for DLAs 
this discrepancy is marginally significant. However, given 
that only 6 of the nearly 600 observed quasar pairs fall in 


quasars is expected to be in halos with masses far from the mean 
halo mass implied from the clustering measurements. 

"" Since we use slices with fixed comoving lengths to mimic the 
velocity windows around galaxies, the width of the resulting veloc¬ 
ity window becomes redshift dependent, but remains close enough 
to the value used in the observational analysis. 
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Figure 10. Cumulative covering fraction of Hi systems with different column densities inside a given impact parameter, R, as a function 
of impact parameter for Lyman-break galaxies. From top-left to the bottom-right, panels show the covering fraction of Hi systems with 
Whi > cm“^, A^hi > 10^’^-^ cm“^, Whi > 10^^ ® cm“^ an d Wht > cm ~^, respectively and with impact parameters 

< R. The data points with error bars show the measurements from [Rudie et ahl ll2012h for a sample of LBGs with M 200 ~ 10^^ Mq 
at (z) = 2.4 and 2.3 for the inner and outer impact parameter bins, respectively. Predicted mean covering fractions for halos with 
IQll-S ^ M 200 < Mq in the Ref-L100N1504 EAGLE simulation are shown with dark-colored regions which indicate the systematic 

uncertainty in the mean due to uncertainties in the background ionizing radiation and the redshift range of the observed LBGs (see the 
text). The light shaded areas indicate the 15-85 percentiles for the scatter due to object-to-object variation. The predictions agree well 
with the observations. 


this bin, one may question the robustness of the error es¬ 
timates. There could also be biases. For example quasars 
covered by DLA absorption may be missing from the bright 
sample because of obscuration. The theoretical prediction is 
also most uncertain at the smallest impact parameters. For 
example, radiation from local stars thought to be the dom¬ 
inant source of_Jiydrogen_£lmtoioima^n close to galaxies 
fe.g.. ISchavell200a : iRahmati et al.ll2013bl) and would reduce 
the abundance of Hi (see Appendix . The presence of 
bright quasars will strengthen this effect. Quantifying the 
impact of local radiation on our results would require de¬ 
tailed radiative transfer simulations that also account for 
the duty cycle of quasars. 

The top panel of Fig.[9]shows the simulation predictions 


from iFumagalli et al.l (l2014l) using open squares and light- 
red shaded areas which, respectively, show the mean and 
1(7 scatter for covering fraction of LLSs around 5 simulated 
galaxies with halo masses M 200 ~ Mq at z = 2. Their 

LLS covering fractions are significantly lower than both our 
predictions and observations. There are several potent ial ex¬ 
planations for this difference. IFumagalli et al.1 ll2014h anal¬ 
ysed the Hi distribution at lower redshift and with lower 
masses than the objects in the QPQ6 samp l e. In addition, 
the simulations analysed bv IFumagalli et ^ ||2014|) did not 
include the efficient feedback that, as we show in lj5] is re¬ 
quired to obtain reasonable stellar masses and Hi covering 
fractions for the halos they consid ered. Furthermo r e, be - 
cause they used zoom simulations, IFumagalli et al.l (|2014|) 
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only considered the distribution of absorbers with small line- 
of-sight separations from the galaxies {R ~ r 2 oo) when cal¬ 
culating covering fractions. In contrast, observations used 
a line-of-sight velocity window of AV = 3000 km/s which 
translates into distances much larger than the virial radii 
of the relevant halos. While this difference does not affect 
the covering fractions at impact parameters R r 2 oo, it 
results in large differences at R > r^oo (e.g., up to > 100% 
difference in the covering fraction of LLSs at i? 200 —1000 
pkpc for M 200 ~ 10 ^^ ® Mq halos at z = 2.5, as shown in 
Fig.lCTI). 


4.2 Hi distribution around LBGs at z • 


In addition to quasars, there are observational constraints 
on the Hi distribution around star-forming galaxies at 2 : ~ 
2 (e. g., iRudie et al.l I 2 OI 2 I : iRakic et al.l 12012 : |Turner_e£_al 
[ 2 OIJ). Here we compare with the data from iRudie et al 
ll2012il . who used spectra of background quasars behind a 
sample of LBGs at z ~ 2 — 2.5 to measure the covering frac¬ 
tion of Hi close to LBGs, which have halo masses M 200 ^ 
10 ^^ Mp (lAdelberger et Si l2005l : iTrainor fc Steidell I 2 OI 2 I : 
iRakic et alT 2013h . By considering all absorbers that are 
within a line-of-sig ht velocity w i ndow of AV = 1400 km/s 
around each galaxy. [Rudie et al.l (l2012l l calculated the aver¬ 
age covering fraction of Hi systems with four different ab¬ 
sorption strengths, and within impact parameters 100 pkpc 
and 200 pkpc for samples of 10 and 25 galaxies respectively. 
Noting that the typical virial radius of those galaxies is ~ 90 
pkpc, the chosen impact parameters are close to one and two 
virial radii. 

To compare the EAG LE predictions with the observa¬ 
tions of iRudie et al. I ll2012ll . we select all simulated galaxies 
with 10^^ ® < M 200 < 10^^'^ M© at z = 2.2 and a = 2.5 
in the Ref-Ll DON1504 simulation, which results in nearly a 
thousand halos at each redshift. This redshift range closely 
matches the redsh ift distribution of the galaxies used by 
IRudie et al.l (I 2 OI 2 I I. We then calculate the Hi covering frac¬ 
tions by adopting velocity windows of AV — 1350 and 
1312 km/s around galaxies at z = 2.5 and 2.2, respectively. 
To account for uncertainties in the strength of the UVB radi¬ 
ation, we adopt the same approach as in the previous section 
and recalculate the Hi distributions afte r reducing the Hi 
phot oionization rate of our fiducial model dHaardt &: Madaul 
l200ll ~) by a factor of 3 (i.e., to Tuvb = 7 x 10“^® s“^ at z 
= 2.5). To bracket the redshift range of observed galaxies 
and the range of possible UVB photoionization rates, we 
take the covering fraction results based on our fiducial UVB 
at z = 2.2 as the lower limit and the lower UVB model 
at 2 = 2.5 as the upper limit. Fig. [10] shows the predicted 
cumulative covering fraction profiles where dark colored re¬ 
gions show the area between the mean covering fractions of 
the two bracketing cases, and the light-colored areas around 
them show the range between the 15 percentiles of our lower 
limit and the 85 percentiles of the upper limit, indicating the 
object-to-object variations in the covering fraction within 
the sample of simulated halos. Counter clockwise from the 
top-right, panels show, respectively, the cumulative cover¬ 
ing fraction profiles of systems with Ahi > 10®® ®cm“^, 
LLSs (Vhi > lO^'^-^ cm"®, sub DLAs (Vm > 10®® cm"®) and 
DLAs (Vhi > 10^®'®cm“^). Note that the quantity shown 
on the vertical axis, /cove(< R), is different from the previ- 
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Figure 11. Differential covering fraction of LLSs around (and 
within AV = 1294 km/s from) halos with mass 10 ^^'® < M 200 < 
1912.2 ^ — 2,2 as a function of impact parameter for 

simulations with different amounts of feedback. The solid blue 
curve shows the L = 25 cMpc reference model, Ref-L025N0376. 
Also shown are simulations using the same box size, resolution, 
and initial conditions but without any feedback {NoFB; black 
dotted), without AGN feedback {NoAGN; purple dot-dashed), 
with half as strong stellar feedback { WeakFB, green short-dashed) 
and with twice as strong stellar feedback (StrongFB, red long- 
dashed) as model Ref. For each model the typical stellar mass 
of the central galaxies in the halos is indicated in the legend. 
While the covering fraction is substantially lower in the absence 
of any feedback, all the models with feedback predict similar Hi 
distributions even though the stellar masses vary greatly. 


ous plots and indicates the total covering fraction of systems 
with impact para meters smaller tha n R. In each panel, the 
measurements of IRudie et al.l (|2012l) at impact parameters 
R = 100 pkpc and R = 200 pkpc are shown as filled circles 
and the error bars show the statistical 1 —cr errors. Note that 
these errors cannot be compared directly with the intrinsic 
1 — a scatter due to object-to-object variation (15-85 per¬ 
centiles) in the predicted values (shown by the light-coloured 
areas). 

As figure shows, the predicted coveri ng fractions 
agree very well with the observed values from IRudie et al.l 
ll2012ll . 


5 IMPACT OF FEEDBACK 

There are two ways in which feedback can affect the com¬ 
parison between the simulated and observed distribution of 
Hi around galaxies. First, feedback can change the distri¬ 
bution of gas around individual galaxies. Second, feedback 
can change the stellar mass of galaxies that reside in halos 
of a fixed virial mass. This will affect comparisons with ob¬ 
servations of the gas around galaxies of a fixed stellar mass, 
since the Hi covering fractions are sensitive to halo mass. 
In this section we investigate the effect of feedback by com¬ 
paring simulations that use different feedback implementa- 
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Figure 12. The same as Fig. llll but for galaxies with stellar mass 
10®'® < M* < Mq. The typical halo masses corresponding 

to each stellar mass are indicated in the legend. The distribution 
of Hi is highly sensitive to the strength of the feedback if the 
galaxies are selected by stellar mass, with more efficient feedback 
yielding higher covering fractions. 

tions, a box size L = 25 cMpc, and our default resolution 
(i.e. N = 2 X 376® particles for this box size; see Table [T]). 

Fig. [11] shows the impact of feedback on the average 
differential covering factor of LLSs around simulated halos 
with 10®®'® < M 2 00 < 10®^'^ M(7) a t 2 ; = 2.2 (i.e., galaxies 
simila r to LBGs : lAdelb ereer et ^ l2005l : iTrainor fc Steid^ 
I 2 OI 2 I : iRakic et al.r 2 OI 3 I 1 . Shown are the reference model 
{Ref-L025N0376, blue solid), a model for which both stel¬ 
lar and AGN feedback are turned off {NOFB, black dotted), 
a simulation without AGN feedback {NoAGN, purple dot- 
dashed), and models with respectively half and twice the 
amount of stellar feedback as Ref ( W eakFB, green short 
dashed; StrongFB, red long dashed). See lCrain et al.l (l2015l l 
for more information on these simulations. The typical stel¬ 
lar mass of the central galaxies in the chosen halos is shown 
next to the name of each simulation. Varying the strength 
of stellar feedback by a factor of two has a dramatic effect 
on the stellar mass, which increases by nearly an order of 
magnitude from StrongFB to WeakFB. However, the effect 
on the Hi distribution is small, with stronger feedback yield¬ 
ing slightly higher covering fractions in the inner haloes. On 
the other hand, turning off feedback altogether does lead to 
a large (up to a factor of 2) reduction in the covering frac- 
ti on of LLSs. The simila rity of our NoFB results to those 
of iFumagalli et ^ (l2014h suggests that the reason why they 
found low Hi covering fractions may be that stellar feedback 
is highly inefficient in their simulations. 

As Fig. [ 11 ] shows, the stellar mass does not change 
significantly between NoAGN and Ref. This indicates that 
AGN feedback does not have a very strong impact on galax¬ 
ies with the halo mass range chosen for this figure (10®®'® < 
M 200 < 10®^'® Mq at 2 = 2.2). However, AGN feedback does 
become important for more massive galaxies, e.g., boosting 
f<r 2 oo Fy ~ 20% for galaxies with M 200 ~ 10®^'® Mq at 
2 = 2.2. 


Because of the sensitivity of stellar mass to feedback, 
the LLS covering fraction is sensitive to feedback if the stel¬ 
lar mass is held fixed, as would be appropriate when com¬ 
paring to observations of the gas around galaxies selected by 
stellar mass. This is shown in Fig. 1121 where we compare the 
mean differential covering fraction of LLSs around galaxies 
with stellar masses in the range 10®'® < M<, < 10®°'^ Mq 
at 2 = 2.2 (i.e., galaxies similar to LBGs; IShaolev et al.l 
I 2 OO 5 II . The curves correspond the models presented in the 
previous figure except that StrongFB is missing because at 
2 = 2.2 this simulation does not contain any galaxy with 
M* > 10®'® Mq. The typical masses of the halos in the 
which galaxies reside are shown next to the model names. It 
is evident that the covering fraction increases rapidly with 
the efficiency of the feedback. Hence, the distribution of Hi 
around galaxies selected by stellar mass, is sensitive probe 
of th eir halo mass (see also iKim &: Croft! I 200 I : iRakic et al.l 

I 2 OI 3 II . 

Comparing Figs. fTT] and 1121 we see that in contrast to 
the other models, the covering fraction predicted by model 
Ref is nearly the same for the samples selected by halo and 
stellar mass. This suggests that (only) this model reproduces 
the stellar mass - halo mas s relation of LBGs, w hich is con¬ 
sistent with the finding of dFurlong et al.l[2014h that model 
Ref agrees with the observed galaxy stellar mass function 
at these redshifts. This means that when we compared the 
simulations results with observations in EH and EH we 
could have chosen to match the stellar masses of the sim¬ 
ulated galaxies to the observed values instead of matching 
their halo masses, without obtaining different results. 

We note that the sensitivity of the distribution of 
Hi to feedback depends somewhat on the column density. 
Stronger absorbers, e.g., DLAs, are slightly more sensitive 
to the feedback effic i ency, consistent w i th previous studies 
dTheuns et al.l I 2 OO 2 I : lAltav et al.l l2013l : iRahmati fc Scha^ 
I 2 OI 4 II . 

We conclude that the inclusion of a relatively efficient 
stellar feedback is necessary to increase the covering fraction 
of LLSs around LBGs to the observed values (see 114.21) . At 
fixed halo mass, but not at fixed stellar mass, the results are 
insensitive to the precise efficiency of stellar feedback. AGN 
feedback on the other hand, has a mass dependent impact 
and helps to boos the covering fraction LLS around bright 
quasars to the observed values (see EH- 


6 SUMMARY AND CONCLUSIONS 

The observed high covering fractions of strong Hi absorbers 
around high-redshift galaxies and quasars has been iden¬ 
tified as_a^_cha]lengefo£simulntions_of_galaxy_f 2 HS®;tio'® 

fe.g.. lFumagalli et al.ll201 j : iFaucher-Gigufere et al.ll2015l ). It 
is therefore important to test whether the EAGLE cosmo¬ 
logical, hydrodynamical simulation, which, thanks to the im¬ 
plemented efficient stellar and AGN feedback, reproduces a 
large number of observed galaxy properties over wide ranges 
of mass and redshift, can also reproduce the Hi observations. 

We combined the EAGLE simulation with photoion¬ 
ization corrections based on radiative transfer simulations 
of the UVB and recombination radiation, to study the dis¬ 
tribution of relatively high Hi column densities (i.e. LLSs; 
Ahi > 10®®^cm“^). Because the main EAGLE simulation 
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uses a 100 cMpc box size, it includes a statistically repre¬ 
sentative sample of galaxies, with a relatively high resolution 
for a simulation of this kind. 

We first demonstrated that EAGLE reproduces the ob¬ 
served column density distribution of strong Hi absorbers 
(i.e, LLSs and DLAs) at 2 = 1 — 5 and yields an evolution 
of the cosmic Hi density that is in agreement with observa¬ 
tions. We then analysed the Hi distribution around galax¬ 
ies from 2 = 4 to 2 = 1, bracketing the era during which 
the cosmic star formation rate peaked. We found that the 
mean covering fraction of LLSs within the virial radius of 
galaxies, f<r 2 oo^ evolves strongly from ~ 70% at 2 = 4 to 
< 10% at 2 = 1. However, the LLS covering fraction depends 
only weakly on halo mass at a fixed redshift, particularly for 
M 200 10^^ Mq. We also showed that /<r 2 oo i® insensitive 

to the specific star formation rate, which suggests that the 
distribution of LLSs is regulated on time scales that are 
longer than the typical time scale for episodic fluctuations 
in the star formation rate. 

At a fixed impact parameter from galaxies, the covering 
fractions of LLSs, sub DLAs and DLAs increase rapidly with 
halo mass and redshift. However, for a fixed redshift and 
after normalising the impact parameters to the halo virial 
radii, the covering fraction profiles of strong Hi absorbers 
around galaxies with M 200 10^^ M© depend only weakly 

on halo mass. The covering fraction profiles of strong Hi 
absorbers in and around massive halos are thus nearly scale- 
invariant and have characteristic lengths similar to the virial 
radius. Exploiting this result, we presented a fitting function 
that reproduces the covering fraction profiles for each class of 
strong Hi absorbers (i.e., LLSs, sub DLAs & DLAs) around 
halos with M 200 ^ 10^^ M© at different redshifts. 

For a given halo mass and redshift, there is a signifi¬ 
cant intrinsic scatter around the mean LLS covering fraction, 
which is related to the complex geometry of the gas distri¬ 
bution around galaxies. This relatively large scatter limits 
the predictive power of studies that use small numbers of 
galaxies to model the distribution of Hi. 

We compared our predictions with measurements of 
the covering fraction profiles of strong Hi abso rbers around 
LBGs and bright quasars at 2^2 — 3 from iRudie et al.l 
(l2012ti and IProchaska et al.l (l2013bl L finding agreement. 
This success may be due to our use of a cosmological sim¬ 
ulation that includes the efficient stellar and AGN feedback 
that is required to produce galaxies with properties close to 
those of the observed population at different epochs. In ad¬ 
dition, instead of choosing a fixed mass, redshift and virial 
radius for calculating the Hi distributions, we found it to 
be important to match the observations more closely. We 
matched not only the redshift range and halo masses, but 
also the line-of-sight velocity interval used in observations 
for finding the absorbers. Moreover, we compared our pre¬ 
dictions with the observed covering fractions at the impact 
parameters that are probed by the observations instead of 
normalising them to the virial radii which is problematic 
since the observed quasars span a range of redshifts and, 
presumably, halo masses. 

Noting that earlier studies showed that local sources of 
ionizing radiation m ainly affect Hi absorbers that are v ery 
close to galaxies le.g.. lSchavel2006| : lRahmati et al.ll2013bll . it 
is likely that they would only significantly change the Hi cov¬ 
ering fractions for R r 2 oo- Therefore, we do not expect the 


neglect of local sources in the present study to change our 
main findings, although it may explain our over-prediction 
of the Hi covering fractions at the smallest observed impact 
parameters from bright quasars (see Fig. O. Accounting for 
the impact of local sources of radiation on the distribution 
of Hi absorbers requires complex modelling of sources in 
addition to accurate radiative transfer. We postpone such 
analysis to future work. 


We also found the assumed strength of the UVB ra- 
diation to be importan t. While our fiducial UVB model, 
iHaardt fc Madaul (l200ll '). produces results that are in rea¬ 
sonable agreement with the observations, we found that 
using a model wi t h a t hree times weaker UVB, similar t o 
iHaardt fc Madai] (l2012ll and iFaucher-Giguere et al.l ll2009ll , 
improves the agreement with the observed column density 
distribution of LLSs and weaker absorbers at 2 « 2.5. 
We note, however, that differences in the UVB intensity 
cannot explain the discrepancy between previous simula¬ 
tions and both EAGLE and the observations, since the 
earlier models used hydrogen phot oionization rates that 
are close to this lower value (e.g., Fumagalli et al. I I20i4 
iFaucher-Giguere et al.ll2015l : ISuresh et al.ll2015l L 


We tested the impact of feedback on our results by com¬ 
paring EAGLE models that do not include stellar and/or 
AGN feedback and models that use a factor of two stronger 
or weaker stellar feedback. The impact of AGN feedback 
on the Hi covering fraction become stronger with increas¬ 
ing halo mass and helps to boost the LLS covering frac¬ 
tion around bright quasars. While efficient stellar feedback 
is required to increase the Hi covering fractions to the ob¬ 
served values, varying its efficiency by a factor of two does 
not change the results significantly at fixed halo mass. This 
suggests that the Hi distribution around galaxies is mainly 
determined by the cosmic supply of neutral hydrogen into 
halos. This conclusion is consistent with the lack of a strong 
correlation between LLS covering fractions and specific star 
formation rates, the rapid evolution in the Hi distribution 
around galaxies, the scale invariance of Hi distribution in 
and around massive halos and the filamentary structure of 
Hi systems around galaxies. 


However, at fixed stellar mass the Hi covering fractions 
are highly sensitive to the efficiency of the feedback. Of 
the EAGLE models analyzed here, only the reference model 
matches the observations of LLSs around LBGs both when 
the galaxies are selected by the halo mass and by the stellar 
mass corresponding to the observed galaxies. This confirms 
that the fiducial EAGLE model reproduces the relation be¬ 
tween stellar mass and halo mass for these galaxies. 


We have shown that a careful comparison with the ob¬ 
served covering fraction of strong Hi absorbers, matching 
the mass and redshift distribution of the observed galaxies 
and quasars as well as the allowed velocity differences be¬ 
tween absorbers and galaxies, results in agreement between 
EAGLE and the data. We conclude that these observations 
therefore do not point to a problem in our general under¬ 
standing of galaxy formation. Noting that EAGLE was not 
calibrated by considering gas properties, its success in repro¬ 
ducing the Hi distribution around galaxies was by no means 
guaranteed. 
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APPENDIX A: IMPACT OF VARYING THE 
UVB 

Self-shielding against the ionizing background radiation 
starts at Ahi > 10^®cm“^. Systems with lower col¬ 

umn densities may not be dominated by neutral hydro¬ 
gen. As a result, the abundance of those systems is not 
only sensitive to the distribution of hydrogen, but also 
to the intensity of the UVB radiation. Observational con¬ 
straints on the intensity of the UVB at 2 < ^ < 6 are 
mode l dependent and uncertain within a factor of a few 
(e.g., iBaitlik et al.l fl98^: iRauch et al.l I l997l: iBolton et alJ 


2003: Faucher- Giguere et al. 20081 : Calverlev et al. 2011 : 
Becker fc BoltorJl2013ll . Models for the UVB are also un¬ 


certain due to the assumptions they need to adopt (e.g., the 
escape fraction of ionizing photon s) and differ from each 
other by a factor of a few (e.g., Hamdt_^_Mada^ 1200 ll : 
iFaucher-Giguere et al.l l2009l : iHaardt fc Madaul 20121) . This 
in turn makes the predicted abundance and distribution of 
Hi absorbers with Ahi < 10^® cm“^ uncertain. 

The impact of varying the UVB photoionization rate 
on the Hi GDDF at z = 2.5 is shown in Fig. lAll The 
dashed red curve shows the resu lt using our fiducial UVB 
model llHaardt fc Madaul 1200 ll ). The long-dashed green 
curve shows the Hi GDDF calculated by assuming a con¬ 
stant UVB photoionization rate at all densities (i.e., no self¬ 
shielding) which, as mentioned above, starts to deviate from 
the reference result at Ahi 10^®cm“^. The solid blue 
curve shows the result of reducing the amplitude of the hy¬ 
drogen photoionization rate by a factor of 3 from that of our 
fiducial UVB model, which improved the agreement with the 
observed abundance o f Hi absorbe r s with column-densities 
Ahi < lO^’^cm”^ from iRudie et al.l (l2013l ). As Fig. [^^illus¬ 
trates, this weaker UVB model produces somewhat higher 
covering fractions for LLSs and weaker absorbers, which is 
also in better agreement with the observations (see Fig. HOD . 


APPENDIX B: IMPACT OF LOCAL STELLAR 
RADIATION 

In I Rahmati et al.l (l2013bl) we post-processed a cosmologi¬ 
cal hydrodynamic simulation with full radiative transfer of 
the ionising background, recombination radiation and lo¬ 
cal stellar radiation. Here we use those simulations to es¬ 
timate the impact of local ionizing radiation from stars on 
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Figure Al. CDDF of neutral gas at 2 : = 2.5 for the EAGLE 
Ref-L100N1504 simulation and different UVB models . The ob¬ 
servational data points are identical to those shown in Fig. ^ 
The dashed red curve shows our predict ion using the fiducial 
UVB model, i.e., iHaardt &: Mad^ ll200lh while the solid blue 
curve shows the result of using a 3 times smaller hydrogen pho¬ 
toionization rate. The long-dashed green curve shows the result 
of using the fiducial UVB model without any self-shielding (i.e., 
the optically-thin limit). The ratios between the two CDDFs and 
that of the fiducial model are shown in the top panel. The obser¬ 
vational measurements, in particular th e grey star-shaped data- 
points at Nm < 10^^ cm“^ taken from iRudie et al] lj2Q13ll with 
(z) 2.4, are in better agreement with the model with a 3 times 

weaker UVB radiation. 



Figure A2. Differential covering fraction of LLSs around galax¬ 
ies with M 200 > 10^^ Mq at 2: = 2.5 using different UVB models. 
The short-dashed orange curve shows the result of using our fidu¬ 
cial UVB model l|Haardt &: Mad^l20Qll L The solid blue curve 
shows the result for a UVB model in which the hydrogen pho¬ 
toionization rate is reduced by a factor of 3 compared to our 
fiduc ial UVB model and is consistent with the jHaardt Sz Mad^ 
l2012h model. The long-dashed green curve shows the result of 
not taking into account self-shielding. For calculating the cover¬ 
ing fractions only absorbers within a line-of-sight velocity window 
of Ay = 3150 km/s around galaxies are taken into account. The 
difference between different UVB models affect the distribution 
of LLS around galaxies significantly and is an important source 
of uncertainty in the predicted covering fraction profiles. 

shown) the reduction due to local sources varies from « 10 
per cent at r 2 oo to « 60 per cent at 0.1 r 2 oo- 

Finally, we note that local AGN could potentially also 
have a large impact. This is, however, even more uncertain 
because the ionizing radiation may be anisotropic and vari¬ 
able in time. 


the distribution of Hi around galaxies. The underlying cos¬ 
mological simula tion uses the OWLS reference sub-grid feed¬ 
back model (see lSchave et alll2010l l in a periodic box with 
L = 6.25 cMpc, using cosmological parameters consistent 
with Wilkinson Microwave Anisotropy Probe 7 year results 
and a resolution s imilar to that of the re ference model in the 
presen t work (see Rahmati et al j|2013bl for more details). 

In lPahmati et all 1 2013bl 'l we showed that photoioniza¬ 
tion by local stellar radiation becomes more important than, 
or comparable with, ionisation by the UVB in regions less 
than r 2 oo away from galaxies. 

To illustrate the impact on the Hi covering fractions, we 
show the differential covering fraction of LLSs around galax¬ 
ies with M 200 ~ 10^^ Mq at « = 2 in Fig. IBII once without 
including the local stellar radiation (i.e., in the presence of 
the UVB and recombination radiation; red solid curve) and 
once after including it (blue dashed curve). The impact of 
local stellar radiation on the covering fraction of LLSs is 
strongest very close to galaxies, the reduction of the cover¬ 
ing factor due to local sources increases from « 10 per cent 
at r2oo to « 20 per cent at i? ~ 0.1 r2oo- For DLAs (not 


APPENDIX C: ALLOWED LOS VELOCITY 
DIFFERENCE BETWEEN ABSORBERS AND 
GALAXIES 


As discussed in 1 12.31 when we associate Hi absorbers with 
galaxies we take into account their relative line-of-sight ve¬ 
locities. The typical velocity differences used in observa¬ 
tional studies are AV ~ ±1000 km/s, which corresponds 
to cosmic scales that are much larger than what has been 


(e.g.. 
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2014: 

Faucher-Gieuere & Kere^ 

I2QI1I; 
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shows how LLS covering fraction predictions change by vary¬ 
ing the size velocity window that is searched for absorbers 
around galaxies with M 200 > 10^^ Mq at 2 ; = 2.5. The 
dashed orange curve corresponds to the value we used to cal¬ 
culate the Hi cover ing fraction profi l es for c omparison with 
the observations of IProchaska et HI ll2ni3bh who sued a ve¬ 
locity window of AU = 3000 km/s around galaxies. While 
increasing the size of the allowed velocity window does not 
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Figure Bl. Differential covering fraction of LLSs around galax¬ 
ies with M 200 ~ 10^^ MfTi at 2 = 2 in the full radiative trans¬ 
fer simulations of iRahmati et aT] ll2013lJl with and without local 
stellar radiation shown by blue dashed and red solid curves, re¬ 
spectively. Local ionizing radiation from young stars reduces the 
covering fraction of LLSs close to galaxies by ^ 20% while at 
larger impact parameters, R > r200i the impact of local stellar 
radiation on the LLS covering fraction becomes smaller. 

affect the result for R r 2 oo, it increases the LLSs cover¬ 
ing fraction for R > r 2 oo- Note that the size of the smallest 
velocity window we showed in this figure, AV = 450 km/s, 
is still « 5 times larger than the velocity window that corre¬ 
sponds to the common choice in previous theoretical studies 
that considered the regions confined within the virial radii 
of halos with M 200 10^^ M© at z = 2.5. 


APPENDIX D: NUMERICAL CONVERGENCE 
TESTS 

To study the impact of numerical resolution on our results, 
first we use two different cosmological simulations that have 
identical box sizes of 25 comoving Mpc, but different resolu¬ 
tions. The first simulation, Ref-L025N0376, has a resolution 
that is identical to that of our fiducial simulation (i.e., Ref- 
LIOONI 5 O 4 ) and the second simulation, Ref-L025N0752, 
has an identical box-size and sub-grid physics but 8 times 
better mass resolution. As Fig. IDll shows, the LLS cover¬ 
ing fraction around halos with 10^^'® < M 200 < 10^^ M© at 
2 ; = 2.5 is not fully converged with resolution and increases 
by increasing the resolution of the simulation. The dotted 
green curves shows the results in the Recal-L025N0752 simu¬ 
lation, for which the feedback implementation is recalibrated 
to reproduce similar galaxy properties to those found in the 
Ref-L025N0376. As the figure shows, recalibrating hardly 
has any impact on the covering fraction of LLSs around 
galaxies (for details of the Recal-L025N0752 model see S15). 
However, as the shaded area around the solid blue curve 
shows, the intrinsic scatter around the mean LLS covering 
fraction is larger than the difference between the results at 
different resolutions. Moreover, other uncertainties, such as 



Figure Cl. The differential covering fraction of LLSs as a func¬ 
tion of normalised impact parameter for different line-of-sight ve¬ 
locity differences between absorbers and galaxies. For impact pa¬ 
rameter R > r200 the covering fraction increases strongly when 
the velocity cut is increased. Note that the velocity window cor¬ 
responding to the virial radii of halos with M 200 > 10^^ Mq at 
this redshift is < 100 km/s 

the amplitude of the UVB radiation, have larger effects on 
the LLS covering fractions than the resolution effect. In fact, 
for any resolution, the UVB model should be recalibrate 
such that the cosmic distribution of Hi absorbers is well re¬ 
produced. As the long-dashed curve in Fig. IDll shows, this 
would reduce the differences in the covering fractions of sim¬ 
ulations that have different resolutions. 

The box-size sensitivity of the LLS covering fraction 
around halos with 10^^ ® < M 200 < 10^^ Mq at 2 = 2.5 
is shown in Fig. ID2I where the Ref-L100N1504 simulation 
is compared with simulations with identical resolution but 
factors of 2 and 4 smaller box sizes, the Ref-L050N0752 and 
Ref-L025N0376 simulations shown with green dotted and 
red dashed curves, respectively. The LLSs covering fraction 
at R > r 2 oo increases with increasing the simulation box-size 
but converges for Lbox it 50 cMpc. This highlights the im¬ 
portance of having a large cosmological box for successfully 
simulating the enhanced distribution of LLSs out to large 
impact parameters. 


APPENDIX E: HYDRODYNAMICS 

As mentioned in ^ in EAGLE we used “Anarchy” (Dalla 
Vecchio in prep) for hydrodynamics instead of using the de¬ 
fault hydrodynamics implementation of G ADGET-3 • “An¬ 
archy” uses the SPH formulation derived bv iHopkinsI (l2013f l 
in addition to modified artificial viscosity switch and time 
step limiters (see Appendix A in S15 for more details). 

The difference caused in the Hi distribution by using 
“Anarchy” instead of the default GADGET-3 hydrodynam¬ 
ics is shown in Fig. lEll for a halo with M 200 = 10^^ ® M© 
at 2 = 2.2 and in the absence of any feedback. The left 
and right columns show the result obtained using “Anar- 
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Figure Dl. Differential covering fraction of LLSs around halos 
with < M 200 < 10^^ Mq at z = 2.5 as a function of 

normalised impact parameter for simulations with different res¬ 
olutions. The dashed magenta curve shows the Ref-L025N0376 
simulation, which is identical to our fiducial simulation (i.e., Ref- 
LIOONI 504 ) except for having a smaller box size of 25 comov¬ 
ing Mpc. The dot-dashed blue curve shows a simulation with 8 
times higher mass resolution and the shaded area around it shows 
the 15 — 85 percentiles for the Ref-L025N0752 simulation. The 
dotted green curve, which is almost identical to the dot-dashed 
curve, shows a high-resolution simulation recalibrate to achieve 
weak convergence, i.e., reproducing galaxies with properties very 
similar to those in the fiducial resolution. The long-dashed purple 
curve shows the result from the Ref-L025N0376 simulation but 
after modifying the UVB such that the global CDDF of LLSs 
becomes identical to that of the Ref-L025N0752 simulation. For 
calculating the covering fractions only absorbers within a line-of- 
sight velocity window of AV = 3150 km/s around galaxies are 
taken into account. While a higher resolution results in larger 
numbers of LLSs and therefore a higher LLS covering fractions, 
the difference is small compared to the intrinsic scatter of the 
covering fraction and other uncertainties like the intensity of the 
UVB radiation. Indeed, as the long-dashed curve shows, if the 
simulations with different resolutions are normalised to have the 
same cosmological distribution of LLSs, the covering fraction of 
LLS around galaxies becomes insensitive to the resolution. 


Figure D2. Differential covering fraction of LLSs around ha¬ 
los with 10^^'^ < M 200 < 10^^ Mq at z = 2.5 as a function 
of normalised impact parameter for simulations with different 
box sizes. The solid blue curve shows the Ref-L100N1504 sim¬ 
ulation. Dotted green and dashed magenta curves show simula¬ 
tions with box sizes of 50 and 25 comoving Mpc, respectively. The 
shaded area around the solid curve shows the 15 — 85 percentiles 
for the Ref-L100N1504 simulation. For calculating the covering 
fractions only absorbers within a line-of-sight velocity window of 
AV = 1294 km/s around galaxies are taken into account. The 
LLS covering fraction within R > r2oo increases with the size of 
the simulation box but the simulations with Lbox ^ 50 cMpc are 
nearly converged. 


chy”and GADGET-3 , respectively. The Hi distribution looks 
smoother in results obtained by “Anarchy” and the resulting 
LLS covering fractions, f<r 2 oo^ slightly higher that those 
in GADGET-3 . This trend can be seen in the differential LLS 
covering fraction profiles of galaxies with similar masses, as 
shown in Fig. IE3I However, the typical differences in the 
covering fractions are smaller than the variations expected 
due to orientations of galaxies, object to object variations 
and the impact of feedback. 

In the presence of stellar and AGN feedback, the Hi 
distribution produced by “Anarchy” and GADGET-3 look 
slightly different as shown in Fig. lE2l and Fig. IE3I The LLS 
covering fractions, however, are almost identical. We con¬ 
clude that the impact of different hydrodynamics implemen¬ 
tations on the Hi covering fractions, which was small in the 
absence of feedback, becomes even smaller in its presence. 
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Figure El. The impact of various hydrodynamics formalism on 
the distribution of Hi around a simulated galaxy in a halo with 
M 200 = Mq at 2 : = 2.2 and in the absence of feedback. 

The left column show the Hi distributions around the halo calcu¬ 
lated using “Anarchy”, our fiducial hydrodynamics implementa¬ 
tion and the right column shows the same but using the standard 
GADGET-3 implementation. The virial radius of the halo is in¬ 
dicated with the blue circle centred on the halo. The size of the 
region is 500 x 500 pkpc. The LLS covering fraction, /<r 2 oo’ 
dicated on the top-right corner of each panel. While the covering 
fraction of LLSs, /<r 2 oo slightly larger in the results obtained 
using “Anarchy”, the difference is smaller than the variations ex¬ 
pected due to orientations of galaxies, object to object variations 
and the impact of feedback. 
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Figure E2. Similar to Fig. lEll but in the presence of stellar 
and AGN feedback. Despite producing slightly different Hi dis¬ 
tributions, both SPH implementations result in very similar LLS 
covering fractions. The impact of different hydrodynamics imple¬ 
mentations on the Hi covering fractions, which was small in the 
absence of feedback, becomes even smaller in its presence. 



Figure E3. Differential covering fraction of LLSs around halos 
with < M 200 < Mq at z = 2.2 as a function of im¬ 

pact parameter for simulations with different hydrodynamics and 
feedback implementations. The solid blue curve shows the REF- 
L025N0376 simulation which uses the fiducial feedback imple¬ 
mentation and “Anarchy” hydrodynamics implementation. Long 
dashed red curve show the result of using standard GADGET- 
3 SPH implementation in the presence of our fiducial feedback. 
The black dotted and purple dot-dashed curves show simulations 
without any feedback which use “Anarchy” and GADGET-3 , re¬ 
spectively. The median stellar mass corresponding to halos in each 
model is indicated on the left-hand-side of the relevant name. For 
calculating the covering fractions only absorbers within a line- 
of-sight velocity window of AV^ = 1294 km/s around galaxies are 
taken into account. The difference in the LLSs distributions cased 
by varying the hydrodynamics is much smaller than the difference 
caused by feedback, and the typical scater due to galaxy to galaxy 
variations. 


© 0000 RAS, MNRAS 000, 000-000 






















